{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Robust Linear Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation\n",
    "\n",
    "Load data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sm.datasets.stackloss.load(as_pandas=False)\n",
    "data.exog = sm.add_constant(data.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huber's T norm with the (default) median absolute deviation scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-41.02649835   0.82938433   0.92606597  -0.12784672]\n",
      "[9.79189854 0.11100521 0.30293016 0.12864961]\n",
      "                    Robust linear Model Regression Results                    \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                   21\n",
      "Model:                            RLM   Df Residuals:                       17\n",
      "Method:                          IRLS   Df Model:                            3\n",
      "Norm:                          HuberT                                         \n",
      "Scale Est.:                       mad                                         \n",
      "Cov Type:                          H1                                         \n",
      "Date:                Tue, 24 Dec 2019                                         \n",
      "Time:                        15:02:54                                         \n",
      "No. Iterations:                    19                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "var_0        -41.0265      9.792     -4.190      0.000     -60.218     -21.835\n",
      "var_1          0.8294      0.111      7.472      0.000       0.612       1.047\n",
      "var_2          0.9261      0.303      3.057      0.002       0.332       1.520\n",
      "var_3         -0.1278      0.129     -0.994      0.320      -0.380       0.124\n",
      "==============================================================================\n",
      "\n",
      "If the model instance has been used for another fit with different fit\n",
      "parameters, then the fit options might not be the correct ones anymore .\n"
     ]
    }
   ],
   "source": [
    "huber_t = sm.RLM(data.endog, data.exog, M=sm.robust.norms.HuberT())\n",
    "hub_results = huber_t.fit()\n",
    "print(hub_results.params)\n",
    "print(hub_results.bse)\n",
    "print(hub_results.summary(yname='y',\n",
    "            xname=['var_%d' % i for i in range(len(hub_results.params))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huber's T norm with 'H2' covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-41.02649835   0.82938433   0.92606597  -0.12784672]\n",
      "[9.08950419 0.11945975 0.32235497 0.11796313]\n"
     ]
    }
   ],
   "source": [
    "hub_results2 = huber_t.fit(cov=\"H2\")\n",
    "print(hub_results2.params)\n",
    "print(hub_results2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Andrew's Wave norm with Huber's Proposal 2 scaling and 'H3' covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters:  [-40.8817957    0.79276138   1.04857556  -0.13360865]\n"
     ]
    }
   ],
   "source": [
    "andrew_mod = sm.RLM(data.endog, data.exog, M=sm.robust.norms.AndrewWave())\n",
    "andrew_results = andrew_mod.fit(scale_est=sm.robust.scale.HuberScale(), cov=\"H3\")\n",
    "print('Parameters: ', andrew_results.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See ``help(sm.RLM.fit)`` for more options and ``module sm.robust.scale`` for scale options\n",
    "\n",
    "## Comparing OLS and RLM\n",
    "\n",
    "Artificial data with outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "x1 = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x1, (x1-5)**2))\n",
    "X = sm.add_constant(X)\n",
    "sig = 0.3   # smaller error variance makes OLS<->RLM contrast bigger\n",
    "beta = [5, 0.5, -0.0]\n",
    "y_true2 = np.dot(X, beta)\n",
    "y2 = y_true2 + sig*1. * np.random.normal(size=nsample)\n",
    "y2[[39,41,43,45,48]] -= 5   # add some outliers (10% of nsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: quadratic function with linear truth\n",
    "\n",
    "Note that the quadratic term in OLS regression will capture outlier effects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.88961659  0.54542326 -0.01414964]\n",
      "[0.45159041 0.06971947 0.0061691 ]\n",
      "[ 4.53587549  4.81389359  5.0871971   5.35578603  5.61966037  5.87882013\n",
      "  6.13326531  6.3829959   6.62801191  6.86831333  7.10390018  7.33477244\n",
      "  7.56093011  7.7823732   7.99910171  8.21111563  8.41841498  8.62099973\n",
      "  8.81886991  9.0120255   9.2004665   9.38419293  9.56320477  9.73750202\n",
      "  9.90708469 10.07195278 10.23210629 10.38754521 10.53826955 10.6842793\n",
      " 10.82557447 10.96215506 11.09402106 11.22117248 11.34360932 11.46133157\n",
      " 11.57433924 11.68263233 11.78621083 11.88507475 11.97922408 12.06865884\n",
      " 12.153379   12.23338459 12.30867559 12.37925201 12.44511384 12.50626109\n",
      " 12.56269376 12.61441184]\n"
     ]
    }
   ],
   "source": [
    "res = sm.OLS(y2, X).fit()\n",
    "print(res.params)\n",
    "print(res.bse)\n",
    "print(res.predict())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate RLM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.81621613e+00  5.30643114e-01 -3.59526003e-03]\n",
      "[0.13815687 0.02132956 0.00188734]\n"
     ]
    }
   ],
   "source": [
    "resrlm = sm.RLM(y2, X).fit()\n",
    "print(resrlm.params)\n",
    "print(resrlm.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare OLS estimates to the robust estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xbdcea58>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsAAAAHTCAYAAAApha5HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XdclfX7x/HXzRKcOHPiFgdu3CMcpZmamqXlaJmmqdmwsmnrZ4pprjT3yspKTb+VmqbmwIGiohm5RXCgCAoc1jn3748rRHIhcjiM6/l4nAd4zg3nc75fjTefc32uyzBNE6WUUkoppfIKJ0cvQCmllFJKqaykAVgppZRSSuUpGoCVUkoppVSeogFYKaWUUkrlKRqAlVJKKaVUnqIBWCmllFJK5SkagJVSSimlVJ6iAVgppZRSSuUpGoCVUkoppVSe4mLvJyhRooRZqVIlez+NUkoppZTK4/bu3XvJNM2Sd7vO7gG4UqVKBAYG2vtplFJKKaVUHmcYxun0XKclEEoppZRSKk/RAKyUUkoppfKUdAVgwzAeMAxj63/uW2MYRgP7LEsppZRSSin7uGsNsGEYRYFFQIEb7usHHDdNc39GnjQpKYmzZ88SHx+fkS/P1tzd3Slfvjyurq6OXopSSimllLqF9ByCswJ9gJ8BDMMoBnwBzDQMo51pmpvu9UnPnj1LoUKFqFSpEoZh3OuXZ1umaXL58mXOnj1L5cqVHb0cpZRSSil1C3ctgTBN86ppmtE33PUq8APwNTDQMIzu//0awzAGG4YRaBhGYERExE3fMz4+nuLFi+eq8AtgGAbFixfPlTvbSimllFK5RUYOwTUEZpimeR5YDvj99wLTNGebpulrmqZvyZK3bsWW28Jvitz6upRSSimlcouM9AE+BlQB/gZ8gXT1W7sfq4LC8F8XQniUhbKeHozu5E2PhuXs/bRKKaWUUioXysgO8ARguGEY24G2wPzMXVJaq4LCGLMimLAoCyYQFmVhzIpgVgWFZcr3Hzt2LJs3b77lY/v372f//gyd81NKKaWUUtlUugOwaZp+/34MN02zi2marUzTfMg0zWt2Wx3gvy4ES5I1zX2WJCv+60Ls+bSABmCllFJKqdzI7qOQ71d4lOWe7k+PK1eu8MQTT2C1WjFNE19fXzp37kxsbCzVqlVjwYIFjBkzhpUrVwKwZMkSNm7cSExMDL17905znVJKKaWUylmy/SS4sp4e93R/esyePZuuXbuyadMmXF1dOXfuHCNGjGDDhg2cOnWKCxcuMG7cON5++23efvttNm7cCHDL65RSSimlVM6S7QPw6E7eeLg6p7nPw9WZ0Z28M/w9T548Sf369QHw9fXF1dWVuXPn0q9fPyIjI7FYbr27nN7rlFJKKaVU9pXtA3CPhuUY16su5Tw9MIBynh6M61X3vrpAeHl5cfjwYUDqfOfNm0fv3r359ttvKVDg+sA7PDw8iIuLA2TIxe2uU0oppZRSOYdhmqZdn8DX19cMDAxMc9+RI0eoVauWXZ/3Ti5dusQTTzyBaZokJSXRqVMnli9fTtGiRbFarfj7+9OqVSsiIyN58sknsVgsjBs3DoBhw4bddN1/Ofr1KaWUUkrlRYZh7DVN0/eu1+XFAGxvuf31KaWUUkplR+kNwNm+BEIppZRSSuUQpgk2m6NXcVcagJVSSiml1L1LSoLkZPl8/Xrw84OiRWHfPocuKz2yfR9gpZRSSinlYPHxEmyDgmD/fvl46BD8+iu0by87v/Hx8NRTkAMaBWgAVkoppZRSqSIiJOAGBUGrVtC6tYTdlIP/xYtDw4YwYgSUKSP3deoktxxCA7BSSimlVF5kmhAbCwULQlwc9OkjoTcsLPWajz+WAFy3LqxeDQ0aQPnyYBiOW3cmyLM1wNOnT8fPzw8PDw/8/Pyujz1WSimllMqV/voLliyB116Ddu2gWDEYMkQe8/CA6Gip4/3iC/jjD4iMhPffl8fz5YNu3aBChRwffiEb7ACPGiWlJJmpQQP48ss7XzN8+HCGDx9OtWrV2Lx5c+YuQCmllFLKUeLjpWRh3z7Z4X31Vbm/b18IDgZ3d6hXT3Z8O3aUxwwD/vzTcWvOYg4PwNmJn58fTZo04eDBg6xbt46xY8fi5+eHn58fCxcuBODJJ59k4MCBXLx4kbp16zJjxgzHLloppZRSeVdsbOqhs0mTYNEi2elN6c5QqVJqAJ45E4oUgZo1wSVvR0CHv/q77dRmpZ07dzJy5Ej8/f1ve83s2bPx8fFh7Nix9OrVi4MHD1KvXr0sXKVSSiml8qSoKNi7V3Z2U27Hj0vpQoEC0pasbFno2lUOqTVqBJUrp379LabX5lUOD8DZiY+PD7169brlYxaLBQ8PD0JCQtixYwebN28mKiqKsLAwDcBKKaWUylwXLqSG3Oefl24LixZJ7ShAxYoScPv3l+AL8NZbclN3pQH4BgULFkzzZzc3NyIiIgBYu3YtPXv2xNvbm6ZNm/Lcc8/xv//9Dy8vL0csVSmllFK5QcrkNGdnOHwY3n5bQm94eOo1jRpJAO7ZE+rUkd3d4sUdt+ZcQAPwHXTv3p1hw4axceNGiv/7F+3FF1/kueeeY8GCBRQuXJhly5Y5eJVKKaWUyhFME86eTS1j2LtXbh98AMOGSaeFEydksESjRtC4sZzsL1xYvt7LS27qvuX5AHzs2LHrn/+3G4SPjw9/3uJE5PLly+29LKWUUkrlZKYJoaEScAsXhg4d4Nq11ADr5AS1a8vwCG9vua9aNdkFVnaX5wOwUkoppVSmGT8etmyBwECZqAbw6KMSgAsXhoULoUYNqF8f8ud36FLzMg3ASimllFL34vx52dkNDJSb1Qq//iqPrV8Ply7J0IjGjeV242H5Z55xzJqzwKqgMPzXhRAeZaGspwejO3nTo2E5Ry/rljQAK6WUUkrdTmSkTOxq317+PGQIzJ4tnxsG1KoFLVtKyYNhwO+/S3lDHrMqKIwxK4KxJFkBCIuyMGZFMEC2DMEagJVSSimlUhw9CmvWwJ49sHu3HEoDOHcOSpeWnd2aNcHXV7ox/KeDVF4MvwD+60KwJFlJivIg/nQJCtUPxZJkxX9diAZgpZRSSqlsITERDh6UoLtnj/TP9faGHTvg9dehQgVo2hQGD4YmTaBoUfm6rl0du+4Msnd5wskQV6J31iQupAyGs438Nc7j7JFEeJQl054jM+XJABwTE8OAAQOIiIigatWqVKhQgZo1a9K/f//r18TGxtK/f38iIyPx8vJi8eLFGIbhwFUrpZRSKkNMUwJvvnxw5Ag8+6yUNSQmyuMlSsCTT0oA7tkTOneGBx5w6JIzk73KE0wTNm2Sc3/n1rfBcEuicNMTFGp8EmcPGc5R1tPj/l+AHTg+AI8aJX8JM1ODBnecsTxt2jSqV6/OypUr6du3L8uXL+eDDz5Ic82SJUto0aIFb775JoMGDSIwMJAmTZpk7jqVUkoplfkiIqR8Ydcu+bh7N7zxBrzzDpQqBe7uMHKk7PA2aSJT1VI2uQoXTu27m0uklCfc6H7KE6xWWLlSgm9goPyu0H9ENLvy7SHROeH6dR6uzozu5H3f67cHxwdgB9i1axeDBg0CoHXr1pQpU+ama8qVK8eiRYvo2bMnc+fOzeolKqWUUio9LBYZKpGUBH5+8tHLC+LjpR7Xxwcef1yCLsgEtS1bHLrkrHa7MoR7LU+Ij4clS8DfX0qlq1WDr7+GgQPB3b0Iq4JqaReIdLvDTq29XLt2jQIFCgCQP39+rl69etM13bp1w2Kx0KtXL9q1a8fkyZNxdnbO6qUqpZRS6r9+/hk2bICdO+Vd5ORkaNFC6nddXWHOHAnBjRvDvz/v87Kynh6E3SLsprc8IToaZs2SyHb+vJz/++EHqRa5MRr1aFgu2wbe/8qTRxULFy5MTEwMILW+hW/xVsfRo0fp3Lkz+/fvJyIigqVLl2b1MpVSSqm87epV2LgRPvtMDqOlWLBAboUKwejREohXrEh9vH9/aNtWw++/RnfyxsM17SZeesoTzp2Ts4FeXvD221C3rvzfsXs39O6dNvzmNI7fAXaAZs2asXnzZrp168bWrVs5cOAAjRs3TnPN3LlzqV27Ns888ww+Pj7Ex8c7aLVKKaVUHmCa8tEwYPFieZ/98OHU+2vVknIHDw+YNw+KFAGXPBlj7lnKrmx6yxP++Uf+51+8WDbXn3gC3nwTGjXKylXbl2Gm/MWyE19fXzMwMDDNfUeOHKFWrVp2fd47SekCceHCBapXr06FChVYtmwZxYoVA+DZZ5+lV69e9OvXD9M0KVKkCN9++y350zmy0NGvTymllMr2YmKk/diOHXLbuVM+envDt99KsWnz5nJr2hQ8PR294lxv92452LZyJbi5wXPPydnBqlUdvbL0Mwxjr2mavne9Li8GYHvL7a9PKaWUuiemCadOSUlCqVLyPvrDD4PNJo/Xri3T1N58E6pXd+hS8xrThHXrJPhu3iy/ZwwbJk0ycmInuPQGYH3vQCmllFKZy2qV/ljbt8O2bbKze+GCvK/+xhtQrx68+66E3mbNUodMqCyTnCwH2caPhwMHoFw5mDhRSq0LFcr49w0IDWDzqc34VfKjRYUWmbfgTKYBWCmllFL3JzoaAgLkVNRDD0FCArRpIy3JqlSR+1q1ko8AJUvCxx87ds15VFycnB/84gs4eVKmOs+fD/36SdnD/dh4YiOPfPMIybZk3F3c2ThwY7YNwRqAlVJKKXXvVq6E33+XXd7gYHkvvW1bCbn588Ovv0KdOnCLXvsq60VGwowZMHUqXLokpdWTJkH37tIuOaNM0yQwPJA5++aw6MAikmwyAS7RmsjmU5s1ACullFIqB7LZZHzw1q1w4gRMmCD3z50r97VoIYMmWrWScoYUHTs6Zr0qjdBQCbpz5kBsLHTpIq3N2rRJHX6XEVHxUXxz8Bvm7JvDgQsH8HDxoGOVjmw8sZFkWzJuzm74VfLLtNeR2TQAK6WUUupmv/0mY762bYPLl+W+cuXgo4+kFdmiRVK7m5ObwdrJqqAwh09EO3xYfldZtkw25/v2lTOG9epl/HuapsmO0B3M3jebHw7/gCXZQsPSDfmqy1c8XfdpirgX0RrgzJaZ/4M+++yzHDhwAHd3d8qXL4+bmxtDhw6ldevW16/x8/OjadOmTJgwgebNm9O5c2fGjh17n69CKaWUymYsFul/9eefsqM7ezZUqgSnT0uK6t5dShvatJF63pRtwxIlHLrs7GpVUBhjVgRjSbICEBZlYcyKYIAsCcHbt8vBtjVr5PeUoUPhtdfk/9KMCAgN4Jejv3At4Rq/n/idI5eOUMitEAPrD+TFRi/SuGzaOQotKrTI1sE3hcMD8Ki1o9h/fv8dr4lOiObghYPYTBtOhhP1HqhHkXxFbnt9g9IN+LLznUcsT5s2jdatW/Pcc8/x66+/MnTo0JuuOXDgADabjUOHDtG5c+f0vSCllFIqJwgOln5Xu3dDYqIEWx8f6dZQqZK0A3jpJUevMsfxXxdyPfymsCRZ8V8XYrcAbLPBL79I8N2+HYoVgw8/hOHDM/57immazNgzg1FrR2E15fXUKVmHed3n8WSdJynoVjATX0HWc3gATo/o+GhspvQKtJk2ouOj7xiA08s0TWJiYnC7zbHH5ORkjh07hpeX130/l1JKKeUQMTHShmzLFmn02r+/bAsWLy7dGkaOhAcflBreG9uR3c/JqDwsPMpyT/ffj6QkmRkyYYJs1nt5wZQp8MILGZ8CfTH2Iov2L2LOvjkcjTx6/X5nw5l+dfvxfMPnM2n1juXwAHy3nVqQ7fcOizuQaE3EzdmNb3p9c9/b6yNGjCAyMpJu3brRvn37W15TtmxZNmzYcNOYZKWUUirbstkkvFqt4OcnE9aSk6VW19c3tclr2bKy+6syVVlPD8JuEXbLenpk2nPExMgZxEmT5JCbj4+MLe7bF1xd7/372Uwbf5z8g9l7Z7Pq71Uk2ZJo49WGvj59mbhj4vX8lZ0Ptd0rhwfg9GhRoQUbB27M1KLqadOmsW3bNvLly8eBAwdueU3Dhg1ZuHAhTz31FNHR0ff9nEoppVSms1hkh/ePP+RWpAisXSuB19sbWreWINyy5f1NOFDpMrqTd5oaYAAPV2dGd/K+7+8dEQHTpsH06XDlipRlz5wpnR0y0tHhfMx5Fu5fyJx9czhx5QTFPIoxvOlwXmz0IrVKykTbR6o9kiMOtd2rHBGAwT5F1UOGDKFNmzY0bNjwlo83atSIMWPG8Nlnn7F9+/ZMfW6llFIqQ6zW1M4Lb7whiSgxEVxcoEkTCbwp5s51zBrzsJQ638zsAnHypAyumD9fft957DHp6NCy5b19n4DQADad2oSHiwfbQrexOmQ1ybZk/Cr58Wm7T+lZqyfuLu5pvianHGq7VzkmANtD0aJFad++PXPnzmXPnj0ULCgF3e+88w4gAbhevXq4ZuT9BKWUUiozWK0QFJS6wxsQAGfOyE6vj4/U8LZvL8FXd3izhR4Ny2XKgbcDB6S+9/vvpaqlf38YPRpq1ZJuE60+T3/IXh2ymt7Le18fVFEkXxFGNRvFoEaD8C5x/7vTOY1hmqZdn8DX19cMDAxMc9+RI0eoVauWXZ/XkXL761NKKWVHpil1vM7OsHo1DBwoo4YBateWsPvuu1C6tGPXqezCNOW84vjxUslSsCAMGQKjRkH58nLNf1utgZRZjOtVN00Itpk21h9ff72210QynxNOjPUby/sPvp+lry0rGIax1zRN37tdl6d3gJVSSqls4dw52LgRNmyQ28SJcqKpRg148kkJvX5+GnpzMZsNfv5Zgu+uXVCyJHz6qXSqu7E5B9y91dq5a+eYHzSfuUFzORV1ipL5S/J03af56chPJFmTcHN2o2OVvD2pz2EB2DRNjPuZwZdN2XtHXSmlVC5gmnJqKSpKShcOH5b7ixeHDh2kQwNAzZoymELlWgkJsHQp+PtDSIjMGvnqK3j2WRlkcSu3aqlmYuP41W30+n46q0NWYzWttK/cnvEdx/OY92Pkc8nHy6Ev58oDbRnhkADs7u7O5cuXKV68eK4KwaZpcvnyZdzd3e9+sVJKqbwjOVlajv3+u+zwVq4sfas8PaFhQylz6NgRGjTQ/rt5xNWrMmn6yy8hPFz+Gnz3HTz+uJxnvJOUVmsJTkeIc9qNzYjB4rQPq9MFtp4pwWstXuPFRi9SvXj1NF+XWw+0ZYRDAnD58uU5e/YsERERjnh6u0oZr6yUUkoBcmpp7lzZ7TUMaNwY6tRJfXzJEsetTWW58+dlWMXMmVLa3aEDLFwov/+kd0/w9YerM3TVp1x0mgvYwIB8tmqMavwhn3V+nnwu+ez5EnKFdAVgwzAeAH40TbPNDff5AJNN03zoXp/U1dWVypUr3+uXKaWUUtlXbKycXlq3Tjo17NghW3lFi8q23sMPS9opXtzRK1UOcOyYlHYvXChd6x5/HN56S2aTpNeFmAss2L+AOfvmcNH5xA2POPFEnd5M7Do0s5eda901ABuGURRYBBS44T4DmARofzCllFJ525Yt8PHHsG2bJBsPDxktfOmSHFr7t7Wmypv27pWDbT/9JL8PPfustG+uXv2uXwpIeeWmU5v4eu/XrDyykiRbEn6V/Him/jN8vu3z61PahrXobtfXkdukZwfYCvQBfr7hvueATUAneyxKKaWUypauXIH16+G332DwYJlEYLPJiK6RI2WXt00b0LMgeZppSlOP8eOl5LtwYRlc8cor6W/kcSnuEgv3L2T23tkcjTx6fUrb4MaDqVmiJgAPVXlID7VlULr7ABuGsdk0TT/DMIoDPyDh93fTNP1uce1gYDCAl5dX49OnT2feipVSSqmsFBsrJ5V++01KG2w2KFYMpk6Ffv0cvTqVjVitstM7fjzs2wdlykj/3iFDZG7J3ew4s4OFBxZyKuoUW05vIdGaSGuv1gxpPITetXvfNKVN3cyefYA/B8aYppl0uw4OpmnOBmaDDMLIwHMopZRSjhEVJd0arFbpxZsvn8yhrVJFyhm6dIGmTVPHEas8Lz5eansnToTjx6V985w5MGCA/PW5myuWK3y05SOm7pp6fVhF71q9Ges3ljql6tzlq1VGZCQAPwhU/zf8NjAM41PTNN/L3GUppZRSWejQIZm6lrLLa7VKyO3bVwo3T5/WMcPqJlFR0rN3yhS4eFH+ykyYAI89dvffj0zTZFfYLmYFzuL7w98Tnxx//TFnw5lGZRpp+LWjew7ApmnWSPn837IIDb9KKaVyFosFtm+X3lMA48bBsmXSjPXtt+GRR6BZs9TrNfyqG4SFweTJ0sc3JgY6d5aODg8+ePdWZtcSrvFN8DfMCpzFgQsHKOhWkGfrP0vz8s0Z+svQ64fa/Cr5ZclryavSXQOcUb6+vmZgYKBdn0MppZS6q7Nn4Zdf4H//kxNKFouM3qpRQ963zp9fijaVuo2//5aJbUuWyJsEffrI4bYGDe7+tUHngpgVOItvgr8hNimW+g/UZ6jvUJ6u+zSF8skvWAGhAXqo7T7ZswZYKaWUyv5sNmlL5u4uwbdrV7m/UiUYNEj+XLGi3Fe1qsOWqbK/nTvlYNvPP0tN7+DB8PrrMtDvdgJCA1h/Yj2JyYlsOLmB3WG78XDxoK9PX4Y0HkLTck1vmoark9qyjgZgpZRSuUdsrLQpW70afv1Vtudef13alX3+OXTrBrVqpX/klsqzTFNKwsePhz//lHkm770HI0ZAyZJ3/tplwct4ZtUzJNuSAahYpCJTOk9hQL0BFPUomgWrV3ejAVgppVTOZ7PJaK21a+VIfpEi0q0h5b3pokWlSFOpu0hOhu+/l8NsBw9ChQpS7ztoEBQsePuvS7QmsuLICmYFzmLL6S3X73cynBjceDAjm43MgtWr9NIArJRSKmcxTSnGXL1aTiNNnQpOTnJQbfBgOYLfpg246rBSlX5xcTBvnnS8O30aateGRYvgqafu/FfpVNQpZu+dzbygeVyMvUhlz8oM9R3Kgv0LSLIm4ebsRrtK7bLuhah00QCslFIqZzhwAJYulULMo0flvmbNZMvOxQUWL3bs+lSOdPkyTJ8O06bJ561ayeePPiq/V92K1Wblt2O/MTNwJr8d/Q3DMOhaoytDfYfycNWHcTKcGFBvgB5oy8Y0ACullMqe4uNljmybNlLSsGGDNFxt107Ga3XvDuXLO3qVKhOtCgrDf10I4VEWynp6MLqTNz0alrPLc50+DZMmwdy5svvbrZtUybRqdevrA0IDWPPPGi7HXWbt8bWciT5DmYJleK/te7zY6EUqFKmQ5no90Ja9aQBWSimVfVy7JofXVqyQjzExsuvbr58UYb74IhQu7OhVKjtYFRTGmBXBWJKsAIRFWRizIhggU0PwoUNS3/vtt/Lnp5+Ws5J1bjNzwjRNZuyZwai1o7CasjbfMr5MengS3b274+qspTY5kQZgpZRSjmWzyXvN4eHSVyoxEUqVkmTSq5fs+ILsAqtcy39dyPXwm8KSZMV/Xch9B2DThG3bpKPDL79AgQIwfDi8+ip4ed36a6Lio1h8YDGzAmdx5NKR6/c7G870qtWLx2s/fl9rUo6lAVgppVTWCw2FlSvlVqGC1O+WLQsffght20KLFnefJatylfAoyz3dnx42G6xZI8E3IABKlICPP4aXX4ZixW79NfvO7WPmnpksO7SMuKQ4mpVrxntt3uOLgC90SlsuogFYKaVU1pk/H2bNgj175M916sjY4RTvvOOYdSmHK+vpQdgtwm5ZT497/l6JifDNNzK17cgRmX0yfTo89xysDwmj2+y0dcadfIqx/PByvgr8it1hu8nvmp+nfZ5maJOhNCrTCIAu1bvoobZcRAOwUkop+zl2DFatkkNrLi7Svgxg3Djo2RO8vR27PpVtjO7knaYGGMDD1ZnRndL/d+TaNZg9W/r2hoVBvXoShJ98Uv76pdQZR1kPEe8STMzV0jy/agZJv20kJimamiVqMrXzVAbUH4Cnu2ea762H2nIXDcBKKaUyV0gI/Pgj/PCDtC4DaN0amjeX4KulDeoWUup8M9IF4uJFaQc9YwZERYGfn3R36NQp7dA//3UhXLEe5ILbu0AyGIDpRHFrW9Y88yEPVnzwpvHEKnfSAKyUUur+JSaCmxvs3Cn1uyDjhydNkgltKSeNNPyqO+jRsNw9HXg7cQImToQFCyAhAXr0kFZmzZrdfO25a+f4K2Y+0W4rwZARxZgGhZMfp1D8M1rXm8doAFZKKZUxf/8tM2OXL4fOnWWEVpMmsg2nPXqVHQUFycG2H36Q36kGDoTRo2+uqDFNkz9P/8lXgV+x4sgKkl2TcbPWINHpBGDDwIX8tqYZqjNWOZsGYKWUUvdm2jSZGXvggLy/3KYNNG4sjzk7w7Bhjl2fypVME/74Q4Lv77/L5OvXX5fy8rJl0157NeEqSw4s4avAr/gr4i+KuhdlZNORVM3fg2nrY4lKPES8UzDutrp4OvvcU52xyh00ACullLqz0FBYu1YGURgG7N0L+fPLVLbevW9OH0plIqtV5qKMHy9/9R54QErJX3oJPP89pxYQGsDmU5spX7g820O3s/TgUmKTYvEt68v87vPp49OH/K75AShbMAz/dW6ER9Wy+7Q5lX0Zpmna9Ql8fX3NwMBAuz6HUkqpTHbunBxk++472LFD7gsJgRo1IDlZjtQrZUfx8dIe2t9fmolUqyZlDgMHgrt76nVbTm3h4aUPk2hNBMDN2Y1+dfsx1HcoTco1cdDqlaMYhrHXNE3fu12n/wVTSimV1u+/y/F504S6deHTT6FPH0kgoOFX2VV0NMycCV9+CRcuSHXN8uUyFPDGM5Rnos/wdeDXTNk15Xr4NTB4s9WbfNLuEwetXuUU+l8xpZTKy+LiZFTWsmXQvj288oq0K3v/fQm9tWs7eoUqjwgPl9A7a5b0833oIeno0L59aisz0zTZeHIjM/bMYHXIakzTpGWFlgSGB5JsS8bN2Y0u1bo49oWoHEEDsFJK5UXr18PSpTKKOCZG6ng7dpTHChWCjz5y7PpUnvHPP1LmsHixVNc88QS8+SY0apR6TVR8FIv2L+KrwK/45/I/lMhfgjcflNEdAAAgAElEQVRbvskQ3yFU8qx0vQZYp7Sp9NIArJRSeYFpwl9/yehhkMQRGAh9+0K/ftLJQXv0qiy0e7ccbFu5EvLlgxdekK4OF90CWHdqMwmhfuR3zc9Xe75iafBS4pLiaF6+OUt6LqF37d64u6QWAuuUNnWvNAArpVRudviwlDcsWwZnzsDZs1CmDMyfD6VKSfJQKouYJqxbJ8F382bp4vDOOzBihHR3CAgNoMPiDiQkJ4ABNtOGu4s7T/s8zctNX6ZRmUZ3fQ6l0kMDsFJK5UaBgTB4sEwMcHKSgsqPPoLCheXxChUcuz6VpyQny0G2CROkfXS5cjI35cUXpeIG4OzVs7z7x7tYki1yhwldqnVhSa8lFPMo5rjFq1xJA7BSSuUGcXHw889QujS0aye7vC4ucqqob1/ZXlMqi8XFyZjiiRPh1CmoWVPefOjXTyZnm6bJppObmbFnBqv+XoXVtOJkOAGQzzkf77V9T8OvsgsNwEoplVPZbLBli5we+uknOTrfr58E4HLlpMhSqUy0KigM/3UhhEdZ7jhEIjJSJmJPnQqXLkGLFvK7WLdu8obEtYRrzN2zhBl7ZvBXxF8U8yjGay1eY6jvUM7HnNcDbcruNAArpVRO1aWLFFQWKiQT2QYOhLZtHb0qlUutCgpjzIpgLElWAMKiLIxZEQxwPQSHhsKkSTBnDsTGwqOPSiuz1q2lldmRiCN8tecrFh1YxLXEazQu05gFjy2gT50+eLh6AFC5aGUNvsruNAArpVROcOkSfPutTGf77TcZRfzii/Dss9C9u/xZKTvyXxdyPfymsCRZ8V8XQnW3ckyYIGctTROeflqmttWtC1vPbGXAyq8JuRRC4LlA3JzdeLLOkwxvMpym5ZpipDT5VSoLaQBWSqnsKikJ1q6FhQtlWEVSEjRsKJ0catSAxx939ApVHhIeZbnpvvizRQn6qSo+Y+R3sGHD4LXXoGJFiIiNYOj/PuDrvV9jYmJg8FLjl/io3UeUKlDKAa9AqVQagJVSKrtJSJD2ZMHBsrtbqpT0iXrmGahXz9GrU3lUWU8PwqIsmCZYjpfi6s6qJIQVwyV/ImPHwvDhULw47AnbwwerpvPdoe+ujygGcDKc8CripeFXZQsagJVSKju4fFlKHBYulJA7f77s9q5fD35+4Orq6BWqPO7V9t6M/Owyl3ZUJulSIZwLx1Gq019M+cCTnk2Ls/zwcqavnM7usN0UdCvIi41epGWFlgxaPYhEayJuzm74VfJz9MtQCtAArJRSjrVxI8yaBatXQ2KihN4W/x4AMgzp36uUA8XEwNy5MGlSOc6FliP/AzEU6bqfas2v8EK7ghyKX8bIybOJiIvAu7g30x6ZxsD6AymcT3pOV/asrF0dVLajAVgppbLasWNQpYr0g1qzRlqZvfyylDjUr+/o1SkFQEQETJsG06fDlSvSYGTWLCjic5CFBxZw7PIxXtqwFROTbjW6MbzpcDpU7nDToTYdU6yyI8M0Tbs+ga+vrxkYGGjX51BKqWwvNlY6OMybB1u3wqZNUtoQFSWnh9zcHL1CpQA4eVKmtM2fDxYL9Oghrcx8GsXw8ZaPmbhjIiaSHfrX7c8n7T+hkmclxy5aqX8ZhrHXNE3fu12nO8BKKWVPV67AmDFS33v1KlSrBuPGQe3a8rinp2PXp9S/DhyQUcXffy9vTvTvL63MXB84xozdM+g8aQHRCdHXr3c2nKldsraG37wqIUE60pw+DWfOpL1FRso49mxMA7BSSmW2y5elzKFZMyhYUIZV9OgBL7wAbdpIba9S2YBpSgXO+PHSca9gQRg1Cka+YuOvhPW8sXsavx39DWcnZ56o/QR+lfwYtXaUHmrLC65elbcDTp1KDbk3ht3z52/+mtKlwcsLKleWMw3Z+J0tDcBKKZUZbDY50DZvHqxcKT8ITp6U7g1Hj4KL/udWZR82G7wz+TIzprgQE1oE14KJ9HvZwmfvwapTC+m4agZHI49SumBpPnjwA4Y0HkKZQmUAqFuqrh5qyw3i4yXcnjol/6367y0yMu31Hh4Sbr28ZMRfyudeXtL4uXx5ad+YQ2gNsFJK3a+ffoI33pAfJMWKyXvHL7ygPXtVtpOQAEuXwoefJhF2yhXnWn/g1vxnPEq6Y3P/i3jXP4i3xtGifAtGNB3B47Ufx805++7iqbuIjJR3o1Jux4/Lx5Mn4dy5tNe6uUGlSrJ7+99bxYpQokSOePdKa4CVUspebDb4/XeoVUt2P/Lnlx8S48ZJqYO7u6NXqFQaV6/C11/D5MmSewqUjafwU4u5WuMlLCRjMQDTmZJmB7a9+H80LtvY0UtW6WGacOHCrUPusWNyyPZG5ctD1arQqdPNIbdMGSn+ziM0ACulVHqFh8OCBdIU9dQp+OAD+OgjeOQRuSmVzZw/D1OmwMyZEB0NHTrA9HlXeGHbG0S7LAcjWS40DQonP06B+IEafrOjq1fhn39SbyEhqZ/HxKRe5+wsu7XVqsHTT0vYrVZNbpUrSxmDAjQAK6XU3Zmm/DD54QewWqF9ezk19Nhjjl6ZUrd07Bj4+8OiRXIWqXdv6DX0EJtipjEgaClxbnG4WquQZMYDNgxcyG9rQllPDUgOk5wMJ06khtsbP9544MwwpFTB2xtat4bq1VNDbsWKOjUynTQAK6XUrYSHw2+/SS2vYUDJkvD66zBokPzAUXnWqqAw/NeFEB5loaynB6M7edOjYTlHLwuAvXvld7Mff5SSzoHPWmnYZw0/nJnKU39uwt3FnX51+1GncB9mb7QSlXiIeKdg3G118XT2YXQnb0e/hNzPYpFQe+RI2ts//0BSUup1JUpIyH3kEflYo4Z8rFJFy6wygR6CU0qpFDYb/PGHvF/888+y23vypOy2KIWE3zErgrEkWa/f5+HqzLhedR0Wgk0TNmyQ4LtxIxQpAs8Ni6RQ23ks/nsGp6NP41XEi2G+wxjUaBDF8xe//lqya5DPFa5ehb/+ktuNQffkSfk/DaTmtmpVOU+QcqtZU37JLlbMsevPodJ7CE4DsFJKAezfD08+KS3LiheH55+HwYPlbUWl/tXq8z8Ii7LcdH85Tw+2v90+S9eSnCwNSMaPh6AgOcPUdsRSjhadyqFLB0i0JuJXyY8RTUfQ3bs7Lk76pq9dWCwSbA8dSnsLDU29Jl8+2cGtXTtt2K1eXXdzM5l2gVBKqTsxTdizB+LiZCRx5crS0eGDD6RgUn8oqVsIv0X4vdP9GXWn3VmLBRYuhIkTpWS0Rk0rw6atZofLJ3x/IQguyJS2xT0WM6D+gExdV56WlCS/IP836B4/Lu8egdSd1KoFbduCj48E3tq15b8vzs6OXb9KQwOwUipviY2VscQzZ8K+fdCqFWzbJu8bb9jg6NWpbK6sp8ctd4Az8/DYf8sswqIsjFkRTMxVgzPbyzJlCly8CI1aXmHgW/PYEjedry6fpki+IhgYmMg7u2evns20NeU5EREyG/rgQfl44ICUMqTU6Do5ye5tvXpyQNbHR27VqunQmxxC/19SSuUdU6fKDm90tPywmjFDhlYolU6jO3nfsgY4Mw+P+a8LSfP9k6/lI3xPFZ6ZUApbIrTu+RfNOk9l46Ul7DsXx4MVH2RSp0mULFCSTks66Zjie5GcLAfSUkJuSuC9cUhEmTISdB9+WD76+Eidrr5LlKNpAFZK5V7JybBmDTz4oBwoKVFCRngOHSo7vzlgqpHKXlLKEOx5eCylnCLpcgGid1Ul9nA5wEq+dktp0OsbtkVswD1CujmMaDqC+qXrX//ajQM36pji24mNlXC7b58UTQcFya5uQoI87uoq5QoPPQT160vYrVcPSpVy7LqVXeghOKVU7hMRIcMqZs6UgyjTp8PLLzt6VSoPu5eOC/WGBXJsQ3kslpNQ41dcy53FVnEdVpdzlC9cnpebvMygRoMokb9EFr+KHCQyMjXkpgTekJDU7gslS0KDBnKrV08Cb82a2kM3F9BDcEqpvCc5Wfr0fvutdP/v2BGmTYOuXR29MpWH3a6mF1J3lE1T2k6PHw/Bf/pCg6Xw1LNgWEkywM2swmu+s/i/R17Qbg7/deECBAamBt19++D06dTHvbygYUN46in52LAhlCun7wDlcen6V2QYxgPAj6ZptjEMwwtYDNiAY8AQ097byEopdTvx8bBrl5Q5uLhI783Bg2HYMDmNrZSD/bemF8CSZMV/XQiP+pTj++9hwgQIPmSjRLP11Px0Cn8nr73haieerP0kEx4dkrULz44uX5ZpH4GBqbeUdmOGIQfTWrSQf/+NGknYLV7csWtW2dJdA7BhGEWBRUCBf+8aAgw1TfOIYRi/AXWBg/ZbolJK3cKZMzBrlpQ6REbKn8uWlcaourOjspFbtUizJTrz98YHqD4LTp+LoXSnxZTpO5VzSSG4uJdmUPVBLA1eSpI1CTdnN4a16O6AlTtYdLTs5t4Ydk+cSH28Rg1o0wZ8feXWoAEUKuS49aocJT07wFagD/AzgGma797wWHHg0n+/wDCMwcBgAC8vr/tfpVJKpTh2DN56C1atkj936wbDh8tJbdDwq7KdG1unWS2uXNtbiWv7KmHLF4Z7z9cpUGUe563RNCnZBP9mS3mizhO4ObvxfMPn7+lAW46e7JaYKB0Ydu2S2+7dUrObonJlCblDhsjHxo2ldaFSGZTuQ3CGYWw2TdPvhj/3ATqbpvncnb5OD8Eppe5bfLwcbKtQQd7ubNIEnn1WujlUrOjo1Sl1R6uCwnhj/lHO/X2ZuKQjEFcU53q/YKuwDifDoHft3rzS7BWal2+OkcFf4LLjiObbMk0ZB5wSdnftktrdlG4MpUtDs2by7zxld1fLGFQ62fUQnGEYVYA3gI4Z+XqllEqX8HDp5PD11/L25vr1EoLPntVm8ypHCA6GnyaV4/j2o9D/KXBKAgPyuRTileZvMazJMMoXLn/fz3OnOmOHB+CoKNnRvXF3NyJCHvPwkN3c4cMl9DZrJv/G9Z0cZWf3/BPk35rgb4HnTdOMzvwlKaXyvKAgmfO6fDlYrdLF4ZVXUh/X8KuyMdOU4YKffw6/bg3HtcVMXAZMItlJpog54cRbrd/ggwc/yLTnzKoRzXdlmvDPP7BjBwQEyMe//pL7DUNajT36aGrY9fHR1mPKITLyU+RtwAuY9u9bNR+aprklU1ellMp7UkaMurrKSOI1a2RXaPhwqFrVsWtTKh1sNvlrO348BJwOJN+DX+L8+nKSSaaVVyv2hO0h2ZaMm7MbD1V5KFOfOytGNN9SbCzs2SNBd8cO2LlTOjUAeHpKR4a+faF5cylp0LpdlU3oIAyllGNdviwlDl99JVtm/ftDTIykicKFHb06pe4qMRG++QYmTEzmb1aS78EpJDywnUJuhXi+4fOMaDqCqsWqEhAaYLcpbVlWAxwaKtvbKYH3wAF5lwak7WCLFtCypdy8vcHJKfOeW6l00EEYSqnsLSQEvvwSFi0Ci0XGj1aqJI8VLOjQpSmVHteuwZw5MHH6Fc6VnYNr1+mQP5RynlUY2exLnmv4HIXzpf4S16JCC7uNJ7bLiGabDY4cga1bJfRu3SrtBgEKFJAShjFjJPQ2by7jxpXKIXQHWCmV9UxTDrWFhMCAATBqFNSp4+hVKZUuFy/Cm9MC+PbADyR6nMLJex025zjaVWrHqOajeLT6ozg7OTt6mfcuMVGGTKQE3u3bpcc2SGeGNm2gdWu51auntfgqW9IdYKVU9pGQAN99J0MrfvlFShsWLpRxpKVKOXp1SqXLiRPgP9FkzsEpWDu8Do1tAHSp3pXPOnxGvQfqOXiF9ygmRsoY/vxTAu+uXdJyEGTIRM+eqaG3ShXtzKByFQ3ASin7uXRJprXNmAHnz8su75kzcvK7YUNHr06pdAkKgs8mxLHi2FJoNgXzob+uP+ZsONOyQsucEX6vXpWgu2WL3PbuheRkcHaWscFDh0rgbdVKfzFVuZ4GYKWUfYSHQ7VqUt/buTO89hp07Ki7SCpHME344w/4aHIYWy1fYTT5GrPmZeqWaETP2h/gv8OfRGsibs5u+FXyc/Ryby0qSsoZUgLvvn1S1+vqCk2bwptvwoMPyoE1rbtXeYwGYKVU5jBNeSv1wAEYORLKloWPPoIuXbS+V+UYViusXAkfzNrDkSJfQqPlGM42ulbtweg2o2jt1RrDMOhcrbPdOjpk2JUrEnQ3b5aPBw7Iv0s3Nzmk9u67EnhbtID8+R29WqUcSg/BKaXuT3Iy/PSTDK4IDJS63mPHwN3d0StTKt02Hw9g8qo/2LYNIov9Cl47cDcK8WLjQbzacgSVi1Z29BJvFhMjJQ1//CG3ffsk8Lq7S8h98EG5NWsmE9eUygP0EJxSyv7+/BOeeQZOnYLq1aXed+BADb8qx4iOhlenr2dBQlcZU9wAijqX5f0OX/JCo7RtzBwuPl4GTaQE3l275BdQV1cJvB9+CO3aSeDNl8/Rq1UqW9MArJS6N+fPSwPU6tXBywsqVIDJk6F7d216r7KFVUFhd+2HGx4OH045xqK/p5JUbxa4/Dum2HDi9bbDeLXFK7f61lkrOVneVUkJvNu3Swh2cgJfX3jjDWjfXg6taUmDUvdEA7BSKn2OHIEvvoAlS6BDB/j1Vxlc8eefjl6ZUtf9dyJaWJSFMSuCARkWERJi8tq0Lfx2ZTJm9TU4NXShZakO7I3cdH1McfvK7R2zeNOEf/6B33+XceCbNknnBpC+uy+9JIG3bVsdKazUfdIArJS6sx07YNw4+N//pLThhRfg1VcdvSqlbsl/XUiaccAAliQr7809zafxG9nr+iWUCcK9eHGer/sO7z38MmUKlbHrmOI7ungRNm5MDb2hoXJ/pUrQp490TmnXDkqWzLo1KZUHaABWSt3MZpPdKGdn2YXauRPGjoVhw/QHscrWwqMs1z83TYg748qV2M1YayyAUucpbtbi7bZf83LrAXi4ph4Ms+eY4jQsFmlNlhJ49++X+z09ZXf3nXdkLLgOnlDKrjQAK6VSJSTA0qXg7y8tzPr0kTHFr72mp8hVpkhPfe79KOvpwdnIeCIv7Ce2wE+Y1Q6BSyLuUW34pucCetbthJGVwdI04dAhWLcO1q6Vrg0JCXJwrVUr+PRTCbyNG8svnEqpLKEBWCkldYazZ8thtvBwaNAAihaVxwoUcOzaVK5xt/rc+xUba1IkOpod0eOgyh4wAZwpnfQmMweOpEe9zAvadxQZKbu7a9dK8A0Pl/vr1JF3UR56SOp49d+WUg6jAVgpBQ8/LC2V2reHBQvkB7S+/aoy2e3qc/3XhdxXAA6/GM/Ls5ax5tJkrMUPYZQsiGkaYJiAScd6Tpm6y3wTqxX27EkNvLt3SxmRp6f8W+rUSW7ly9tvDUqpe6IBWKm86PhxmD4dPvlERqB+9hkULgxNmjh6ZSoXu7E+Nz33w51LJvb/E8HQ+TPZaZsBBS5SyKMew2st4OEmlemy7JHro4qHteie+S/m3DkJvGvXSj3vlSvyS2PTpvD++xJ4mzQBF/0xq1R2pP8ylcpLgoLg88/hxx/lB/Ojj8op8w4dHL0ylQeU9fQg7BZht6znrevLb1cysevwMX4MWsqx/EvAI4GysV0Y2/o1BrVvf72+d+PAjZnb1cFqlZ3dX36RFoBBQXJ/mTLQo4cE3o4doXjx+38upZTdaQBWKi+IjYXevWW3qnBhGD0aXnlFfngrlUVGd/JOE2gBPFydGd3J+5bXp5RMJDgdweIUjDXKk7jEYD532QQe7vgkP8uXvUfRoX7Nm772Xro63HaX+dIlKWn49Vf5txMZKQfVWraU1oCPPCL9ebVcSKkcRwOwUrmVacrwitq15bBNgQLyQ3voUG2irxwipXQhvV0gwqMsWIxgLrq+D0YylATiiuF+bBQHxr1LjfIl7ntNN+4yG6YNz7+DOb12EZGRhyh2KEj+HZUsCV27QpcuUi+fckBUKZVjaQBWKrexWqXEYdw4+PtvOHUKSpeW+5RysB4Ny6XrQFpYVARXr/5GZNGF4JQsd5oGhZ0exadOt0wJvwDT1uynzeEAOhzbTbsTgZSKvQLAX+W9KfbBB1Im1LixjvlWKpfRAKxUbpGYKGOKx4+Ho0fB2xu+/lprElWOsvf03wz/ZjI7LYvhgXiMc00xHwgCw4qBC0VdGt62ZCLdQkNlsuGaNfy6fgP5rElcdcvPn5UbsalqE7ZUacTlAkU5OfbRzHlRSqlsRwOwUrnF6dPw4ovQsKHs9vbooY31VY5gmiYr92/i7Z8ncdT4BZLzUfrCQMZ2fpVS3Qrz4doVhMYFUiG/Lx8/0uveW5rZbLB3L6xZI7eU6WtVq7KqeXdWVWjEnvJ1SHZO/ZFY7jYH85RSuYMGYKVyqqgoaWV26hTMnQvVq8vJdD2Uo3KIRGsi0zd9z7jNk7jkuh/iSlLz6li+eGooXR4sdf26no1G3Ps3j4uTYRRr1shu7/nzUsbQsqW8S9KtG9SsSb794exfEUxyOg/mKaVyBw3ASuU0ly7BlCkwdapMcOvaFZKSZLRq/fqOXp1Sd7Xu2Do++n0ye8P3kuhyCaJq0cppDjOG9Kd+HfeMf+OICAm8P/8M69dDfLx0PencOfUQ239Kgu71YJ5SKnfQAKxUTrJ2rbQzi4uDxx+Hd9+VscVK5QDHLh/nhe/e4s+In8AAnJx4MOELlrw9igrlM3jI7PhxWLVKQu/27VLu4OUl5UCPPQZt2oCb2x2/RXoP5imlcg8NwEpld6Gh0n+0fn3w9YUnnpA+vrVrO3plSqXL1lM7GL3yC3ZFrwTTkPALODsbdOqYcG/h1zSlnvfnnyX4Hjok99evLxPYHntMfinUMiCl1B1oAFYquzp+XKa2LVokI1W3b4cSJWDBAkevTKm7stqsLA9eyXu/fMGJpJ1gKYrnsbfp27o5i2L7Xh9T7FfJ7+7fLCkJtmxJ3ek9e1bqedu2hS+/hO7doXJlu78mpVTuoQFYqewmJAQ++wyWLZNxxYMHy46vUvfothPO7CgmMYavAuYzfsuXRJonIbIK5c9O47MnnqXfZwVxdoaBoekYU2yxSB3vihWwerUc+vTwkHreTz+V/rwlMqcXsFIq79EArFR2s3Ej/PSTjCp+/XUoW9bRK1I50I0TzgDCoiyMWREMkKkhOCA0gM2nNlOnVB1+/zuAOUGzSDCi4ExL6sZMxH/QYzzc0TlNRcJtxxRfuyZjh3/6ST7GxsrUtcceg5494aGHIH/+TFu7Uirv0gCslKPt3w8ffwydOsGQIfD881LnW7Kko1emcjD/dSHXw28KS5IV/3UhmRaAA0IDaLeoHQnWBDCR+t6/e9He/XUmjGxB48bp+CaRkbLDu2KF7PgmJECpUtC/vxz09POTDidKKZWJNAAr5Sh790rwXb0aihSB9u3lfnd3uSl1H8KjLPd0/70wTZPfT/zO09+9JOEXAIOKl0ay4bMvqVbtLt/g/Hmp5/3pJ9i0ScZ3e3nB0KHQq5f06tUhLkopO9IArJQjvPUWTJgAnp7w0UcwcqR8rlQmKevpQdgtwm7Z+5hwlmhNZFnwt3z8+xecjAuG2OLg7gqGDcNwgWIVOXQtjGrcYof5/HkJvMuXw9at0s2hRg14800JvY0ba+cGpVSW0QCsVFbZuVOmtRUvLru9hQvDiBHyUalMNrqTd5oaYMj4hLOo+Chm7v4a/61TuZIcDhd8cA6aQ0GPJuRr9AeJ7vtxt9UFW420JRbnzkno/eGH1NBbuzZ8+KGUN9Spo6FXKeUQhmmadn0CX19fMzAw0K7PoVS2tn27lDqsXy+7vR984OgVqTzifrtAnIo6xcRtXzJ33zwSzBg43pHyZ97gw/4P8+nfa8HFdtPXlIqJZHfN6LSht04dqWt/4gntX62UsivDMPaapul7t+t0B1gpe9m6VQLvxo1yoG3CBKlxVCqLZGTCWUBoAEsPLuXQ+RC2nt2EaXWCQ0/hc+01Ph7WgMcekxa8Cz7Pd73EomRMJI+EbOfRkO00OXs4NfR++KGGXqVUtqQBWCl7mThRplR98YV0dyhQwNErUuq2bKaNSQGTeOv3t7Bhk64OwU/jZx3P2FfL07Zt2mqFd5uWYPfEOXQ+tJmmoYdxwuRoyYqEDHmNWiNfgFq1HPZalFLqbjQAK5VZduyQUodp06TWd+ZMOdimfUtVNhafHM83B7/h/7Z8wYmrRyT4GgDOjOjrw9QnyqdeHBUl3Ru++44uGzbQxWrlVEkvprR6ij1NO/LkgIftPmhDKaUygwZgpe7Xzp3yVu/69VLqcOyYBGAdYKGysUhLJLMCZ/HFtqlEJl6Acw1wPjYWo/V4TCMRNxc3nmruJ8Mo1qyB776D336DxEQZO/zmm9C3L5Xq1uVVPcimlMphNAArlVGmKdOpfv5ZRrJOmADDhmmpg8rWTl45yeSdXzIncB7xtlg42plCwW8wqkd7qj55ni92FuNczC4eifCgwqDPYdsGiIuTX+hefhn69oUmTbR7g1IqR9MArNS9OnJE6hsNA+rVgxYtJBgULOjolSl1W4HhgUzYNpEfj/wANmfMg0/zwInXGfN8XV6YBRuPnGb15KV8ELyJh/8JoHBiHJH5C3Pi0d5UGf4CtG4tp9+UUioX0ACsVHrt3Qtjx8L//gdbtkDbtlLzq1QG3G+LsvRIOdg2fddXnL56EiOxMObuN6hxZSTvvVKOvn1MXA/sgfe+ofW8xTwWc4WrbvlZV6Mlq2u3ZUfF+pQuVpDtbdtm6rqUUsrRNAArdTcHDkjv3tWroWhR+OwzaNDA0atSOdiqoLA0QyrCoiyMWREMkCkhOCE5gWXByxi76WPOXDslB9tsrtQO/onPX+lIl+pHcfp2LPgsg6NHwc2NPZUas6q2H5uqNiHBxe3698qM0clKKZXdaABW6k4SEuChh+Tgz8cfwyuv6OQ2dd/814WkmdAGYEmypp2ilgFR8VF8HeXS4agAACAASURBVPg1k3ZM4aLlHFwrDQWcwMmGs7ONfi2/pOvH78CePVLC4+cnY7kff5xPZu3L9NHJSimVXWkAVuq//vkH5syBzz+HfPlklKuPj+z+KpUJbrermtHd1tDoUKbsmsKsPbOJTb4GJzriFLCIbm0M1hV6lCRbIm5WK35zf4GSDaVHdd++UC41bGfm6GSllMruNAArleLkSdnlXbJEgu/TT0PDhtCmjaNXpnKZsp4embLbevDCQSbumMiy4G+x2UzM4D4U2DeKSQ0v06/OQgp8u5KA4olsblQUv3qP0WL9m7cdUJGy82zvumSllMoONAArde0avPEGzJ8Pzs4wYgS8/TY88ICjV6ZyqYzutgaEBrDp1CYK5SvE//75hfXH1+GUXABb4Ms0/6sTk8tvoKmlO04rz8s7Fs88Q4sBA2jRokW62pZlZHSyUkrlRBqAVd6VmAhubjKpbedOGDwY3nknzdvCStlDRnZbt57ZSsfFHUm0JgLglFCUQtveYuTR/Ixy+pES4VPgnCs8+igMGCAf8+XLktejlFI5jWGapl2fwNfX1wwMDLTrcyh1TyIiYPx4+OEHOHxY+vcmJYGrq6NXptRNYhNjWbB/Ae//8T5RCVFyp81g0I4qfL3xJE6mDZo1g4EDoU8fKF7csQtWSikHMgxjr2mavne7TneAVd4RHQ1ffAGTJ8tkq3795GPBghp+VbYTERvBjD0zmLZrOpHxlykYURmXolcxDRtuNpPnomNxevcd2e2tUcPRy1VK5UJZ0a/cUdIVgA3DeAD40TTNNoZhuAIrgGLAPNM059tzgUplirAwqFsXrlyB3r3lsNttDgMp5UjHI48zKWAS84Lmk2CNp2aINzO3ufFk6Em2Vc/P1kdr49f+OVp8/JJOZlNK2Y29+5U72l0DsGEYRYFFQIF/7xoB7DVNc6xhGL8ahvGDaZrX7LlIpTIkPh5275aJbeXKybjiXr2ks4NS2UxgeCD+O/z58fCPGDaDh/aXYtKOc3hf+gdLi/Ywzp/WPXvSOn9+Ry9VKZUH2KtfeXaRnh1gK9AH+PnfP/sBb//7+Z+AL7Dpxi8wDGMwMBjAy8srM9apVPolJcHChbLLe/EinDkjHR0++cTRK1MqjR1ndjBn3xyCLwSz9/xe8ie68spuZ97YlYSbc37yDf4Ep2EDKFCxoqOXqpTKYzK7X3l2c9cAbJrmVQAjtYVOASDs388jgZt6RZmmORuYDXIILjMWqtRdWa3w3XcwdiwcOyYHgxYu1HZmKttJsibx2Z+f8fGfH2Niggkv74Z3/nDjUt0+FP7+WQp2bp2u1mVKKWUPmdWvPLvKyCG4GMADiAYK/vtnpRzv5Ek5Ce/jA6tXQ9euGiBUthKbGMu83bP4YvPnnLFeAhMwwMkGCTWe4P/Zu+/oqKoujMO/m0lCQu9gIFQpUlQEFVA6CohSBBHBXtAPUUE6IoogHUSUqkizIggqiIBUpUiPIFKlhl5DCan3+2OLiFJDkkkm77PWLGAySc6sFSZ7zt3n3dknjSMkR4arfh0RkaTm69MhE1IArwbuBaYAtwHLE3VFItdj3jxYsAB694abb4alS+HOO3U4SFKUw2cO8+GMN/lw4wSO+Z3j3l3QYnNWBtc4TVyAS7p0gTzbvh1BKn5FJIXw9emQ15wD7DjOQtd1qzuOUxD4AfgJqAxUdF037nKfpxxgSRIrVtjQinnzIDQUfvsNsmb19qpELrJjdxiDv3yNT079TKR/PA9tcqiwpCaRuTtzX79aBBX/lUW7FlK9UHUqhVby9nJFRFK9RM8Bdl23+l9/7nIc5z5sF7jHlYpfkUS3Zw+8+ipMnw45c1qm70svQVCQt1cmYlyXtTPHMmDhu0zOsBOPC/XXZyPj0pcIurc9D32b4x9BJJWoXECFr4hIckvQIAzXdfcBkxN5LSKXFxcHHo+NLV692hIe2raFTJm8vTIRAJau/Y6xM95h/cH1rMwVTaZAqLXsNo6t7sNNzerRfplD0aLeXqWIiIAmwUlKd+CA9feuWQNLltiY1+3bNblNktQ1Tz+KiyNu9iz6fduJN/P+gesH5ITiK2pyctUk7nohhFc/hdy5k3FNIiJyVSqAJWU6cQIGDID334eoKHj+eYiMtB1gFb+ShK5p+tHevZwbO5qJS0YwsOQxtoVgiQ4Aroeb763Nl1+GJNoFCl+fyCQiKVxMDJw5A6dPg+va2RuARYvg4EG7//Rpe0xICDz1lHfXew1UAEvKs2YN1KplRfBjj1m7w803e3tVkkZcbvrRkB820mjvGk6OHcHIY7MZerfLwXsg98Gi+P30FG71vjj+0aQLDKT749UTtTvH1ycyiUgyiI2Fo0fh5EkoXtzumzMHNm6E48cv3DJlghEj7OMNG9pjzp278HXKl4fz4Qavv26/s/+pWjUVwCLXLDbWWhtKlIDSpaFJE3jlFbjtNm+vTNKYf085Cok4xKNhc6m2/Uc6lz3OyLscTgW6ZN5VBaa+TZG8NXj2kWPMOJOVvedWERpQgYNHCkBo0q3paveLSBpx+jTs22dTT8/fjh2zlCSwFsIvvrD7jx613dvs2e3vAGPGwNSplpmfJYulKZX4R85v7dr278yZIWNGu91004WPf/qpfc0MGS58PDAw+Z7/DVABLN7lujBlCnTvbpdOtm2zRIePP/b2yiSNCskazIFjpym9fxq5I+ZSa3s484pAx1oOsX5++G9qAou6UqV8OTpPgCMZwuk2bT2RMUXIQhEiIkj09gRfn8gkIpexa5fttoaHW6EbHm63776zYrNnTxg06L+f164dBAdbQXvLLbYrmzv3hdt5o0bBRx9Zgevx/PfrvPLKldd3yy039vy8SAWweM9PP0GXLpbqULo0jBwJ6dJ5e1WSloWH8/HOmWxf8DFNm0UQ5YGx5YC4AJy1z+Es70DzekXpOAfKlrVPuadf0rcn+PpEJpE0x3Vt13XbNpg50wrdnTth927Yu9d6a0uUsMjPtm3tcwICrL82Xz44dcoK4EcftSul/yxuc+a8sAvbpo3dLidnziR/qimVCmDxjl9+gfvug4IFYcIEaNny0u8+RZJafLwNVBk5Eve7b9lbMJ42TYKIOv/qGO/gWdGBl0v24fUP7Uf2n5KjPcHXJzKJ+JzzBW54OHz1lRW354vcnTut7aB2bRvi1Lat7dYWKgQFCsDtt18oYJs1s93bfPksBenfU04rVLCbXDcVwJJ8Nm+GsDD7D33PPfDZZ9brq11f8YYjR2D8eBg9mrjt2/jmrkz065aLNZ6DeKIyWPa0E4+/XyDfD36IumUu/WWSqz2hUbl8KnhFUprTp+H7720n95+3/v3h6aetbaF9eztYVqiQ3apWhTx57PPr1rX+3Jw5rWD+t5tuurjnVhKNCmBJeuHh1qf0ySf2n7xhQyt6W7Tw9sokrXFdWLrU+t6+/pqo2CgmNr2ZAU/nYVvsQQJP5YUFvcl3/AmatFlDltsXcn+xK48pVnuCiA+Lj7eNm02b7Ha+wG3e3Ppsz5y58Lssf35LLGrQAIoUsftuv90OpWXNeukCN316u0myUwEsSefECXsXPHSo7aa1bm2H3bTjK8nt1CmYNMkK3/XriciZiVGvlue97Fs5ELWNgP3lYf5wivs3oktnD82aQUBAJeDqY4rVniCSyrkuHD5sBe4ff9ifRYta76zrQsWKEB1tBWzBglbkZs9un5s7N6xfb48PvsRVn4AAyJYteZ+PXBPHdd2rP+oGVKhQwV11Pi9O0pZ16ywv8HyW7/l3xCLJZeNG/uw5gF0rv2B5vmiyxdzEijolmZpuDRExJ/HfXZvYBV2oElqTrl0c6ta99CaNiPgA17XpouvXw9mz0KiR3V+unP2+Oi99ejuXMmaM/fvHH60Ht1gxSymSFM1xnNWu6161MVo7wJJ44uIsE3DjRtv5vf122LHDmvpFkktsLHz7LQwfDgsWsKeAhzpPucT4AewH9uP81hR+7syDd1Wg8+e2wSMiPuTMGcumBbsKOX26Fb7Hjtl9RYpcKICffdZaHW65BUqWtFaGfx42q1s3edcuyUIFsNw414UffrBIsw0b4K67bHxxunQqfiX5HDxoeZajR1uMUIECdG/8MINLryLGbzc4gOvA8rbkOtqBRfNCKFnS24sWkRu2c6clC61fb7cNG2zQw6lTVsju3WujfJs2hTJlLMOwzD9OtV4t61Z8kgpguTFbtsALL8DixdYD9dVX9iLz76gWkaRw/lDb8OE2UCUmBu6/nyX929DXXczMbd9AbDqI9wfHBTeAnCUKkT54LSVLhnh79SJyPU6ftgJ37VprWejf3/prJ0yAt9+26LCSJS1loUwZ24gJDr70oIjrNH1tuPr8fYwKYEmYuDjL7c2Uyd5dDx9uhXBAgLdXJmnB2bPw+ef2c7duHWTJgtv6f8xqVIY+f05kydYu+EfngCU9cda/QIbKC/ErPYf0fqVJF1xUE9REUrpDh6x4zZTJcrpbt4atW+1NL1jh27q1/fnccxapWaJEkvwOmr42/KKkl/ATkYk+7VGSnwpguT779lmk2dat9qJ00022C6whFpIcdu6EESNsVPbx43DrrcSOGsHX5YPou2Io6xcNw/9Mflg8lNATz3N/k2gWVFxJFJmBphCviDKRFCcyEn7+2Ub+rloFK1faxspnn1nEWK5cUKqU/b1cOTtfEhp64cRq/vx2SyIDZyf9tEdJfiqA5dpERMDAgTBkiF1mfukli4VJl07FryQt17WxoMOG2eE2x2FZy2r8VKsQp2/KweTfB7Fz5p94jpeEReMo42lB106BNGkCHk8Gpq8trUuXIilFRASsXm2FbqlSUL++9evWqWMfL17cWhjKl4c777T7br0Vpk3z2pKTY9qjJD8VwHJ1q1dDvXqWk/joo/Duu9bvK5KUzrc5DBtmfX85ckDnzvzUsCwPzHmKmB3zYQc4R2+Bud9QPX9Dugz0o1ati6PMNEFNxEvOjwN2XXjmGVi+3CaCntemjRXA+fLBggW2s5s1q/fWexnJNe1RkpcKYLk017V2h3z5LBqmRg3o0OHCO3KRpLJ7t7U5fPSRRRbddht88glHGtTm/bAxDJzzLDHxMfbYeD9Kxz7O+E8bU768d5ctkqa5rv3fXb78wi17dpg504rgAwfsgNrjj9vvkfLlbTIo2MerV/fq8q9E0x59kwpg+a9Fi6BTJzhyxKbipE9v6Q4iScV1rQdw2LALlzobN4ZXX2XPrYUYvHwIo0e24VzcWdhZFfL/iuOJJV1AIGO61qB8qHeXL5LmnDkDv/9usZcAzZpZEgvY4bUKFS58DGyYRCqVXNMelTSRvFQAywW//25ZvjNm2M5vr17q75WkFRUFX3xhQfVhYbZj1LEjtG7Nlgzn6PdLfyYNm0RcfDxuWEsyhnXmlealqNxkGetPLaR6oepUCr36uGIRuUEHDlibwpIldvvtN7v/5EnImBGaN7dd3EqVLGfXxxKBkrqVSkkTyU8FsJhff4XKle2FrG9feO21S881F0kMhw7ByJHW6nDokGV2fvQRtGjBmpOb6LO4Pd9smgpx6XBXvUju7R3o+EJBWk2AzJkBKvEgKnxFkkRcnPXd//ILPPII5MljVwHbtrXfERUrwhtv2J/nC90mTby75lROSRPJTwVwWhYRYRmqVataT1a/fnZQ4XxflkhiW78e3nvP4o2io+0ATNu2LC0ezPiwCaz99GNW7f8VJzoz7q9dKHrkNbq9moeWLS1wRESSyP799ib0l1+sf/fUKbs/JAQeftgOQFetaru7/iodEpuSJpKfforTopgYGxfbsyfExlreYoYMdulZJLHFx9uo7KFDLTs6fXoLrn/tNdzixRm8bDCdx3UmnnhwgVUvcceJfnTvkIUGDTRUUCTRnThhPfeLF8O990LDhnDunE1TK1vWDqrdey/cc8+FcfZ589pNkoSSJpKfCuC0xHXhm2+ga1cbZFGjhmX7Zsjg7ZWJj/jnIY4i6WFIZBi3fTPBhqXky2dXGV54gbisWZiycQrvfPgoG4+FWeHrAHh4oVkBRj+R5aIoMxG5QfHxluSzYIH127uujQ7OmNEK4EKFbLhMlizeXmmapKSJ5KcCOC0JC4OmTS18fMYMeOABVGVIYjl/iCPL0YN0XvM9j637kSxRZzhe+jayff45NG1KtJ/LpLBJ9FrQn12nt8KRkrDxTTz3DgInmsCAQJ6pUV0/liI34sAB291dtMguoXzwgf25ZIkdNH37bWtnuPvuC2c9HEfFrxclV9KEXKAC2Ndt2WIvhM8/byHjc+bYzq96uCSRTR83gz7zvuTBTT/j57r8WLwSYys04mDpcsxpUpGPV4+gz8JBHIraC/vuIPDXqbSq2oj2Y/3Y76nHwp1KdRC5Ie+9Z+1t54dNZMhgGx3nLV+uTY8UTEN7kpfjum6SfoMKFSq4q1atStLvIZdw6JD1+I4eDZkywc6dencvic91Ld9z0CCYP5/TgcF8dev9jKvQgL1Z8hDPaU55ZhKffhYRsUdgZ1UyruvG6w3v55U2js5biiRERIRtbMyfb728ixZZb33fvnaIrXp1qFYN7rhDmx2S5jiOs9p13QpXe5z+Z/ias2dtF6B/f/t7q1bw1lsqfiVxRUVZksOQIZYfHRLC8LovMLpELQ6n38NZv9nEcoRIZwWu/2nY+AC5t3Sl2+P38vyHF7edK/xd5BotWgSdO8OqVRZVli6dHVQ7dMh6eLt29fYKRVINFcC+5tAheOcdqFfPDhyVLOntFYkvOXoURo2CDz+0PsPbboOJE+HRR8n3+2GOfTOaA553gXh7/J81ybupLwNeuYvm4/+bjZ+Q8HcVzOLz4uJg9WprWZs3D9q3hwcftHeOHo8VujVr2tCJoCBvr1YkVVIB7AvmzLkQM1WokPX9Fizo7VWJL9m+3a4sjBtnVxbq1rVfyrVqgeOw9ehWvt7aj32e8UC8JTrEe6heuhLzJ9x12bbD6w1/17Qk8WkREfDss9bacPy43VeunMVVgo0XXrLEe+sT8SFK2EzNfvsN6tSx23ffwZEjdr+KX0ksv/5qySHFisGYMdCsmf3czZoFtWuz/tAGGkx8jBIflOTz3z7H3dgYPzcIPzwEBwbSp0X9K565ud7w9ysVzCKpyokTFkv50kvw5pt2X6ZMlsveuLGNCD90CNasgUaNvLtWER+kHeDU6OhR6NTJduOyZoXBg+HllzUqSxKH61qBO2CA9RxmzQpdukCbNjYVCvh17690nvkuiw58D1EZcVZ15NHQdvR4Ow8nMy275kSH6w1/17QkSfU++MD651eutGzejBmhRQv7mONYUoOIJDkVwKlRQIC1PbRrZ/PYs2f39orEF8TE2K7TwIGwYQPkz2+H3J5/HjJlwnVdFvy5gE4z3mX18XlwNjsBa3rywq2v0PWTbOTPf/4LVbrmKLPrDX/XtCRJVfbsgdmzLZnhk08si3fTJvuze3e47z7L4v13c7yIJDkVwKlBbCyMHQuTJ9uLaebMNslNhx8kMZw6BR9/bD2+e/ZAmTJ2sK15cwgIYOnupYyaO5pFW9awO2oDnMpLht8G0a7Ki7T7KuMNvf+63vB3TUuSFO+PP+z1+scfLSEF7M1keDiEhtoBUmXxinidCuCUzHXtcFvHjvaiWqWK9fnmzaviV67bv9MTulfITr2FU2D4cOtHrFbNEh7q1QPHIS4+jt4LetNzcQ9cXHAdMv/egXeq96JVr6C/B0jdqOsJf9e0JElxtm+3YrdGDZuyuWOHtTlUrQrPPGMHRkuVulD0qvgVSRFUAKdUhw7BY4/ZaeBixWDaNJvXrhdPSYB/picUPL6PVrO/oWb3ebjxsTiNG1tP+d13AxATF8PYlZ/TY05fDrubwQUc8HP86NQmO69V8+6bL01LEq+KjbXX5ZkzbYNi2za7f/BgK3Rr17ZzGhkzenedInJFKoBTmuhoCAy0vt64OBg2zE4Jq0dMbsDA2Zspsmcz/1s+hQc2LyHG42FqmVp8V+sxvhz8FABRsVG8v3gcfRb356SzEw7cRu4Tb3O4eF9cYsDxxxNb2rtPRMQbwsNh/36LIYuNtZSG+Hjb9X3tNdvlvflme2xgoN1EJEVTAZxSRETY9LZPP7WYqSxZYMEC7fjKjXFdWLyYvqM7UHXnWiIC0zOqYhPGlW/I4YzZcIAz0WfoM2cMQ1cM4qxnH4TfTeXYD2hY807GbVmLX0xvzvmtJyi+LBMWBlIye7h2YMW3xcVZSsPMmTBjBqxbB7ffDmvXWvvZwoVQurSNHxaRVEkFsLfFxtoBpLfesraHFi1szCyo+JWEi4+3X959+8KyZZTJmJX+1Z7i03IPcCqdzSGO5wwx7kJy9HqaKP8jOHtq8EDGiQzpUJMSJRzu6TefyJg40nEL6eJvASAy/vJDKkRStVOnLIcX4OmnbTPCzw8qV7b/R/XrX3jsnXd6ZYmpkSY3SkqlAtibjh+3Oe7nD7jNmKEXVrkxsbHw1Vc2BnvDBpsMOHw4v5S7n/E/bONE3AbOOiuJOnuaqPS/QGAE/n8+QMuQNxjYpzI33XThSylzV3zen3/C99/bIKHFi+3foaE2je2BB2zIkGImE0yTGyUlUwHsDYcOQe7ckC0bVK8O775rk3604ysJde6cDUYZONBOoZcuDZMmwaOPQkAADYBlx36h38ouQBxkAU94NVoXGUKvYXeQJct/v6Qyd8VnrVplCQ0bNti/S5WCDh3A47F/16jhvbX5kOsddS6SnFQAJ6e9ey38/KuvYONGKFwYRozw9qokNYuIgJEjLcP34EGoWBGGDoUHH7TLt8CWQ7t5ccJAFp4eCU4cOODg4a0WdXizxh2X/dLK3BWfcPYszJtnu7w1alibWf78kDOnDXp56KELB9gkUSXXVSS1WUhCqABODqdP2wG3wYPtcEXbtrb7K5JQR4/C++9bSsjJk3D//dC1q2X5/nUlYe3O7bw4qR8rYyYAkOXwA5y9aS7xxBDoCaT2zdWv+C2UuSupluvaFZFvv4W5cyEy0vp7zxe6efPaIWNJUslxFUltFpJQKoCT2tmzcMsttvvbvLkdpihUyNurktTqwAF7IzVyJJw5Aw8/DN26Qfnyfz9k0e9/0PrLPmz0+xziAwg9/CL9G3Skeb0CLN+7jIU7F1K9UPVrGleszF1JNXbssLSGxo3tTeCoUXZV5LnnoEEDe3OoeLJklRxXkdRmIQmlAjipnI/NSZ/eesvuvtsuT4skxO7d1t/70UcQE2NDUrp2tV7fv0xfvo52095lZ/BUiA+mxMl2DHusPfdXunCyrVJopWsqfEVSPNe1yMhp0+z222+QLh0cO2avu7Nm2QE2na3wmuS4iqTDupJQKoAT24YNVvDOng1Ll0KlShaULnIFl+1h277dEh0mWBsDTz4JXbr8fSl32Z5lDJo7iYV//MaxjEvALzPlz3Zj1NNtqVAqpxefkUgSiIuzwtff3/re27e3Aveee2DQIDtMfD6bN0cO765VgKS/iqTDupJQ110AO46TDfgMyA2sdl33xURfVWp04AD06AFjx0LmzPZifMflDxiJnHepHraPR8+g/O7ZhM7+1qYAvvgidOwIBQoAVgO8MmEEw3e+Ak48ZIDb455ncuuBFAvN6s2nI5K4oqPhp5/gm2/sINtHH9lY+AYNrK+3QQPIk8fbqxQv0WFdSaiE7AA/AXzmuu5njuN87jhOBdd1VyX2wlKVmBjL7z1wANq0sUJYuw9yjf7Zw1b64HZeXvoVD2xZytnAIHj9dbv9FdAbE+Py1oR5DAvrxZmci//+Gh4/D81qFVHxK77j1Clo3dpyek+etGK3fn07wAZ2FUTpDWmeDutKQiWkAD4KlHEcJysQCuxJ3CWlEvHx9sL80EO2Q/fhh3bYrXhxb69MUpl9JyK5bd9mXln6JbW3ryQiXQber9yc8RUasHbgYwCcOePSYfQsxm3vRVTu5finD+G+TK/xS+QYouOiCfQEUr1Qde8+EZEbcfq09e1GRNjBtYwZraWscWNo2hRq17YeX5F/0WFdSYiEFMC/APWBV4E/gGP/foDjOK2AVgAF/rpk61MWL7ZdudWrLWanQQO7JCdyvZYv5/NpPam0ZSXHgzIxsMoTTCz/IKfSZSBf1mAOH4nn1eHfMeVQb2JzryZdpoK8FDqSIY8/Q3BgOpbtefS6Uh1EUpSTJ20C5tSpVvyeOwdlytgkNseBNWt0iE1EkoTjuu71fYLjfAK0dV03wnGc14HTruuOudzjK1So4K5a5SMdElu3QufOduI4f36LNGvR4u+BAyLnXTWYfckS6NkT5s4lKms2PizXiE9uq8eZdHaAx3M2iGxR4awJGISbaz3pzxXl1XLdeKfJEwR4Arz0rEQSwYkTkCWLFbatW1ukX0iIRfo1bQr33nthIpuIyHVyHGe167oVrva4hOwAZwPKOo6zHLgb+CkBXyP1iY+3nd69e6F3b2jX7sJpY5F/uGIw+6ntVvjOnw+5csGAAaT73/8ouvUknlnfcOzwOuIOZOVszh8gZBOZo0vS5c5JdKzbHH8/hbZIKnXihB1gmzwZ5syBn3+2aMhXX4UnnrC/ayNBRJJRQn6j9gXGAQWBZcAXibqilCQ6Gj7+GJ5+2ordiRMhNPTCIQyRS/hPMLvrcvu2tRRo0gV2hNmJ9cGDLdkhQwYAwveuZUNUe9wsMZAVssQWpW/1ybSq8jAeP+2GSSoVHg4vvWSxkDExULCgxULmzm0fL1nSu+sTkTTrugtg13VXAKWv+sDUzHWtzaFzZ9i2zWLNHn/ckh5EruLvAHbX5Z5dYby65Avu3vs7BzNmh6FDoVUrCA4mPh6mfxtFp8/Hs71wFwiOAcAPPzrd/yz/q/LIFb/PVdssRJJbRIQdDvb3h0cftTScnTttp7dZM3sNVU+viKQAuqb6b6tW2QG3n3+2KVuzZkHdO9JmSgAAIABJREFUut5elaQiIVmCKBS2nHa/fEaF8D/YnzEHb973Er9UaciC1+oSHQ0Txp2jx/SPOVC0P5TaS16/0hxzzhIXH0egJ5AahWpc8Xtcsc1CRbAkp9OnreidPNleL6OioFYtK4CDgmD9em+vUETkP1QA/1v79rBpk82Rf+4528kQuVYLFvDt5G7kXLeCfZly0v3+1kwuex+e4CB61LyV/kPO0nfuaE6WHgB3HKBEUBWGNh5HnWK1WL53+TUnOvynzQKIjIlj4OzNKoAl6cXEWPwjwDPPwJQpllX94ou201tJiSQikrJddwrE9UrxKRCnTsHAgdanFhICO3bYZbvMmb29MklNFi2Ct96yP0NCCHuiNa9lLM+u03Hk8s9MvsNFmP3nZ0SWGwQZDnNb5poMbdSD6oWrJejbFe4yk0v9z3WAHf3q39BTEbmk2Fg7vPnFFzB9OqxdC4UKwcqVEBlp6Q06yCYiXpaUKRC+IS4OPvkE3nwTDh60w20vvACFC3t7ZZKa/PyzFb4LFtjhyGHD4IUXuC0oiFdWLKPflB9Zs+EAK0tMgXuPcXeOOgxu8Cb3FLjnhr5tSNZgws/3Gv/rfpFEFR4OffrA11/D4cO2OdC4sb2Ggs5GiEiqlDYL4LlzrdVh/Xq45x6L57nrLm+vSlKTpUut8P3pJ0t1eO89u/wbHExYGLT9YDYLQx6E9LFQDsrlvIdRjYZwV77E+TnrWKfERT3AAMEBHjrWKZEoX1/SMNe13d3oaKhYEQIDYdIkOwvx2GNQr5719oqIpGJpswCeNMkObnz9NTRpolPJcu2WL7fCd84cy/EdNAj+9z/c4PQsXAi9Bh9lQeR7UHkgeGIB8DgeHrm1fqIVv3DhoJtSICTRbNsGn38On30GW7bYQbaffrKf88OHNYZYRHxK2iiADx+Gt9+G55+HcuXg/fct11cv6HKtVq+2dplZsyBnTujfH15+mbigDEyfDr3fO8y6oCE4d38Igaepkq8GKw4sJTY+lkBPINULVU/0JTUql08FrySOF16wzHPHgWrVoEMH2xw4T6+VIuJjfLsAPnfOejLffRfOnIFSpawAzpbN2yuT1GL9etvxnTYNsme38ddt2hAVkJGJE6Hfhwf5M88gnBojcPwjaVryUXrUeIMyucuwbM+ya051EEk2p0/bIbbJk+HTT62nt1YtKF7cWhzy5/f2CkVEkpzvpkBMm2Z5vjt3woMPwoABcMstyb8OSZ22bLGrBl9+CZkyWc9427acdDMzahQMHrOfw8UG4tw5CvyjaFGmBd2rvUHJnJpsJSlQTIydffjsMyt+z56FAgXsdfKOO7y9OhGRRKMUiHXrbGdj7lyoXdvbq5HUYscOeOcdG3sdFARdukCHDuyPyk6nfsuYHPYd0Rm24ff49/j5xfLEbY/TrUo3iuco7u2Vi1zMdW23N1Mm6++tX9+ufj3xBLRsaQeAFVsmImmU7+4AR0XZEAuPJ/m/t6Q+e/daq8zHH9vPTOvW0KULW07kZuBAGLdsOnEPPwJ+seDAg8UfZGidoRTNXtTbKxe52K5d1towcSKUL28H28AOtFWtaqkOIiI+SjvAOrQhlzF9bfjf6QmlPZEM3Tmbm6dOgvh4Owz0xhus3JeP/q1h6rxd+FXth9t0DDjxgKU6VM5fWcWvpCxTp8IHH9gwFrDDbA88cOHjuhImIvI33y2ARS5h+tpwun6znsCIE3RcMZWnV39PutgYdj30CAXe68ucrYXp/wQsWLuTgFp98HttPB4P3JH7flYcmI/rxuLijxNd2ttPRdK62FiYN88OsPn720S2/fuhVy94/HGb0iYiIpfkuy0QIpdQu+dM6s77kla/fkPG6Ei+v6UqQyu3YN+ZyqTfVIZ1u/4kfZ0+nCs5AX+PH8+Xe55y2Z5m8KyjnIjbwDm/9QTFlyWrpwx9Hy6rGDJJXq4LYWHW3vD55zbFctYsG1Jx7pxd+VKuuYikYWqBEPmnqCgYPZovB7xFzrMnmHvz3Qyo9CRrDlQlYkoRYv32k/mhZ/FrNJE4jz+t73iJzvd2Jn/m/NzTbz6RMXGk4xbSxVuSSGR8HANnb1YBLMln7147yPbbbxAQYOk2TzwBNWrYxzWdTUTkmqkAFt8WG2sHgt5+G3btYneR23iuQg/mnbmJs3u34h7ciN8DfaD4ZKIDAmhzx8t0vrczIZlC/v4S+05EXvJLX+5+kUQRHQ0//AAREfDkk3DTTVCwoI3cfvRRyJHD2ysUEUm1VACLb3Jd+OYb6N4dNm2CChU43Ocjev5wN98tCYOW94P/OQDi8eehos8wutE73JTppv98qZCswYRfotgNyRqc5E9D0qCwMBg/3t64HTliw3ueeMLSSb77zturExHxCQqBlCQ1fW049/SbT+EuM7mn33ymrw1P2m/oupb9fNdd0LQpAHuHTuGZ0ivI9/R9zJi7nwyPtLHi1wEcaFLsJb574qNLFr8AHeuUIDjg4ji94AAPHeuUSNrnImlPt25w++0wYoS1NsycCStWqK9XRCSRaQdYksz5xIXImDgAwk9E0vWb9QBJ0zu7fDl07QoLF0KBAmztNo5Ovz3O9Lb+pMu3mcLte7Et+Ati/Pzxd/1xXZdATyDtq7a44pc9v9bz0WkhWYPpWKeE+n/lxsTG2gG2ceOgZ08oWxYaNYJ8+aB5c7U4iIgkIaVASJK5p9/8S7YO5MsazJIuNRPvG/3+O7zxBnz7LW7u3Gxs3J1X1rdiwdJ0ZCm6iQKP9+Z3vy8I8g/i5TtfpkPlDmw/tp2FOxdSvVB1KoVWSry1iFzN5s3wySeW5HDgAOTJA2PGQIMG3l6ZiEiqpxQI8bokPzy2Zw+89RbuhAmcCQhmWPGXGXq8O4dH5yVvmU2Uf7cXa2K+YHtAMO3vbE+Hyh3InSE3ALkz5FbhK8nHda2NITLSprOdO2cpDs8+C/XqWaqDiIgkGxXAkmSS7PDYsWPQty988AFx8S5jQp/graO9OLwlFP9iKwl67FkOZv+RCILpWLkj7Su3/7vwFUk2rmttOWPH2q7v4sUQHAyTJ8Mdd0DevN5eoYhImqUCWJJMxzolLuoBhhs8PHb2LAwbBv364UZEEHbrkzTZ1J0/d92Mp+oYPJUGEhu0jTiCCPE8ytrXhpErQ65EejYi1+jwYWtvGDsW/vgD0qeHZs1s1zc4+OLxxCIi4hUqgCXJJNrhsdhYOyj09tuwbx8bCj3IM+f6sCqsLOnKLyaw5mNEpz/fZ+4hZ3QXAs5VUPEryScuDmJibBjFrFnQoQNUqgQffWSZvZkyeXuFIiLyDyqAJUk1Kpcv4WkJrgvTplk01ObNbMlZiVZ+X7I0vAr1n95E7oot+GHPl+D+FVHm2OfE+P1JSOYqifYcRC5rzx470DZ2LLz2GrRvb/F7FSpAqVLeXp2IiFyGcoAlZVq0CLdSJWjShN17/WjIdMpHLqHw67mo/0lLvs1XikUHvqPxzf8jX9ybOASC64eDP1mc25XRK0nr++/tEFuhQhZhdsstcOut9rH06VX8ioikcNoBlms2fW140mfh/v47bqfOOD/M5HBgPrryMT8EP8WjXbZT7+YnmLjlC4J2B9Gxckc6VO5Argy5mL42nB6zsrHn7CpC01fgnXoPK6NXEt/hw5Drr7aa99+HjRstd/q556BwYe+uTURErotygOWa/HuoBdiBtr4Pl02cYnPfPuLe6IEzYRynnUy8G9+V7wq+SrPXd7Plpl5M/sNyfNvc2ebvwlckyUVHw7ffWi/vwoWwcyeEhMD+/VYM+2sPQUQkJVEOsCSqgbM3X1T8AkTGxDFw9uYbK4AjIojqPRC/oYNxY2IZzqtMrFKHkCZzCc3WmN475hJ0Ooj2lS7O8RVJUgcPwpAhdvjy8GEIDbVhK+nS2cdvuvTYbBERSR1UAMs1SfShFjExRAz+CL9eb5Px7GG+oDkzK71LwTarWL/tQdYdj4Pj0KJMC96r+54KX0l6MTFw5IgVt9HR1uZQrx60agX33w8ej7dXKCIiiUQFsFyTRBtq4bocGDkNp1sX8pzcykKqMavGQCr2yI7/iXfoGzYRF2vL8TgeyuQuo+JXktaePdbiMHYslCkDs2fbju+BA5A1q7dXJyIiSUApEHJNOtYpQXDAxTtg1zvUYvO4pWzOdS95X27CsZP+fHD/98StGsvh10byyOISfPX7VzQr3Ywg/yA8jodATyDVC1VP5Gci8pdFi6BBA0ty6N0bbr8dXnnlwsdV/IqI+CztAMs1SehQC9eFpRO2Et+pC1UOf8MBJy9T64whdEgNwrb14/UfGuNxPLS5qw2d7+nMTZluYtmeZSzcuZDqhapTKbRScjw9SSvO7+oGBcGKFXbr2hWef94KYRERSROUAiFJIi4Ovh9/lHNde9Lk8EiinCDW1u5E1g+a8sGW9xi3bhwex0Or8q3ocm8XQjKFAMkUtSZpi+vCggUwciRMn24H2x5/3EZr+/tDYKC3VygiIolEKRDiFefOwaSPozj69ge8dLQ3mTjF5AcfYtXT+dkes5ofJvfCcRxeLP8iXe/tSr7MF4rbf0ethZ+IpOs36wFUBMv1i421onfECNi0CbJnt2ltlSvbx9On9+76RETEa1QAS6I4cQJGjnDZMeBrupzsQhF2cKBcPeb3qceTK14ndkMsAI1LNub9uu8TmiX0P18jyaLWJG05eBDy5LHUho8+gsyZYcIEaNbMWh9ERCTNUwEsN2TfPnjvPVg7YhnvnG1PV5ZxukhZwt//gv6epYxY0Y4414paj+PhzpA7L1n8QhJErUnaER0N33wDw4dDWBiEh0OmTHbQLVs2b69ORERSGKVASIJs3mznhqoX3MGdgx7lp7OVuTPnDg6MHsybH9bg5rBnGLFyBA8Ue+CaUx0uF6l23VFrknYcOADdu1ts2WOP2YS2t94Cx7GPq/gVEZFL0A6wXJdff4X+/WHBtBP08LzLSHcYniAPhzu1Z8BdMQxf152oA1E8eduTdK/SnaLZi15zqkPHOiUuOW75eqLWJA1wXTh92nZ4DxyAvn2hfn1o3doGVvjpfb2IiFyZCmC5KteFH3+0wnfJohjaBY9mUvDbpD93jGPPPMqgh3LwwcZRRK6OpGXZlrxZ9U2K5Sj29+dXCq10TXFmCY1akzTi1CmYOBE+/BAqVIBJkyy7d88eCAnx9upERCQVUQEslxUbC5Mnw4ABEBbm8kSOWXQp25q1WXbxS+7b+LlJI97f9RVnws7QvExzelTrQcmcJW/oezYql08Fr1xs2zYreseNg4gIKF/eRhSfp+JXRESukwpg+Y+zZ+GTT2DwYNi5Ex4q8jszy7Rn94nZ1GrocM4fXMJgaxjNSjejR9UelM5d2tvLFl8SH299vI4Do0bZ4bZHHrFJbRUrXujxFRERSQA1y8nfjh6Fd96BggWtzrgl1xH+rPcy3+66jUyHltPjxeJE+rucH53y6l2v8lXTr1T8SuKJiIBhw6BkSZg71+7r1Al274bPP4dKlVT8iojIDdMOsLBnDwwZYpGpZ85AoweiGVxkOEUm9eTM+lMMaHc3A3Js4ljUFvwcPxwcAj2BNC/T3NtLF1+xaZO1OUyYYAfcKla8MKEtd27vrk1ERHyOCmAfcr1jhIdNOUjvPvEcDsuDA1SrF8n4uvMpOKw9Z+duZUjLEvQr4eFw1DLqhdajZ/WexMbHXlOig8g1i4uD++6DQ4egeXO7/FDhqlMsRUREEkwFsI+4njHCv/wC7bpFsurnPDgBsWQqt4u7iv5ErxUjyDN7HR/Uy02fp7NxIGYztfPVpmf1nlQOrfz356vwlRty5owlOEydCj/8AAEB8OWXUKyYdntFRCRZqAD2EVcbIxwfDzNmWJTZ0qXgn95Dlns3U6hkGB1XjePhqbMZeWcgD7TIyInAQ1QNqcqX1d+hWqFqXnpG4nP27LHDbGPGwPHjluawb581nd9zj7dXJyIiaYgKYB9xuXHB4UfPMWGCRZlt3Gi1xgcfwNBtP/LU+u+5c8EkRpePonXdQCKCokgXV5ifnviQmoVr4uiwkSSWsDAreF0XHn4Y2raFypV1oE1ERLxCBbCPCMkaTPg/iuD4aA+nwwpwZnURnj4Jt94Kn34KzZpBwLwfeaDXy4wvvZd6jwMO4MaQLfp5SmZsTq0itbz2PMRHREfDlClw8iT873/2A9i7t40rLljQ26sTEZE0LsExaI7jjHAc56HEXIwkXMc6JQgO8BB3JpATi4sTPrImx+eXotjNDj/8AOvWQcs7t+BpXJ+vOtSjbouD9Lqou8HB3xNLp7o3NshC0rhjx2w0ceHC0LKlpTq4ru30dumi4ldERFKEBBXAjuNUAfK6rvt9Iq9HEujWLPkotPle9o2qycllN5O16An6TzjEhlXpqHdPBHTqwLQmpbi90I80fwQCby5Gy5Kd8XPSgeuHnxPA69UaawqbJNyECRAaCt26QalSMHOmNZyrzUFERFKY626BcBwnAPgI+MFxnIau6357ice0AloBFChQ4IYXKZe3bp0dbJs8GTyejDz9JHTsCCVL5ob4eNyxY5k1qj1v3nGSNU2heJaifF6rF81KN8Pj5+HlPQ0VayYJ47qwZAnkzQs33wylS1uPzeuvQ9my3l6diIjIZTmu6179Uf/8BMd5DqgPtAZeAQ64rvvB5R5foUIFd9WqVTe0SLmY68LChVb4zp4NmTLBiy/auaJ8f23gukuWMK/XM7yZfyvLQ6FwcAhv3d+Hlre2xN9Prd9yA2JjLcJs8GBYuRJeeglGjvT2qkRERHAcZ7XrulcNk09IJVQOGOO67gHHcT4F3gUuWwBL4omLg+nTrfBduRLy5IE+feyMUdas9phlq79l3KftWRm5nXWVINQ/B2Pq9uHp258hwBPg3Scgqd/IkdCvn40mLlYMRoyAJ5/09qpERESuS0IK4G1Akb/+XgHYlXjLkUuJioKJE2HQINiyBYoWhVGj4KmnICjowoM+HvgYL8ZMIz4rkAVeL9+GPnUHkc4/nTeXL6ldeDiEhFgv7x9/QKFClqX34IPgl+BztCIiIl6TkAJ4LPCJ4zjNgQCgaeIuSc6LiLBCd+hQ2L8f7rgDvvoKmjQBj+fC49Z+PYwec7oyI/9ZizQDPH4ecmYJUfErCbduHQwcaD90P/0E1atb20OAriSIiEjqdt0FsOu6p4BHkmAt8pf9++H99+1qc0QE1K5tO8C1al18oH7j6lm8NfE5pmTfT9bcfryUpz4Tjs4nOi6aQE8g1QtV99pzkFTKdWHuXCt8f/rJGszbtrV2B1DxKyIiPkGnoVKQrVut7pgwwc4ZNW0KnTrZAK1/2rb3N3qOaclnfhvImAF6BNSmXbvPyJo5N0/uWaZUB0m4qCh44gm7xNCvn52uPN9gLiIi4iNUAKcAq1bZwbapUyEwEJ55Bjp0sGSpf9p1fCe9JjzL+OMLCIyHjhGl6fjal+QsXObvx1QKraTCV67d6dPw8cd2uvKnn6ypfM4cKFkS0ql9RkREfJMKYC85f6W5f3+YPx+yZLFBWa++arGqAMv+2s0tnbs0c1Z8wZhtX+HEu7TZlZMuT48lb80G3n0SknodOGAH2UaMgBMnoEoVOHTIDrvddpu3VyciIpKkVAAns/MRqv37w9q1Vm8MHAitWkHmzBcet2zPMmpOrElUbBSu6+KJh+c3puONGj0IfbvzxafgRK5HWBjcdRfExMDDD9vklLvv9vaqREREko0K4GQSGQnjx1uU2Z9/QokSMHYstGz53yvNxyOP021eN87FngMs2KFd5G0MHD0PcuRI9rWLD1i92n7wHnnEprR16mT5vecPt4mIiKQhKoCT2PHjdpV5wOA4Io57CLzpOCUe30Oftjl4uHy+ix57Ovo07y9/n0G/9OdEzCk88XZ/oH86Hn51pIpfuT6uCwsWQN++1t9buLDt+Ho80KuXt1cnIiLiNSqAk8jevfDeezBmjJ0zylD0GHnqbiNd6DHOOfDG9H34+UGjcvmIjIlk5KqR9P25D0cij9JgM7yzNhtn277MwluCqV64hg62yfVZutTiy1autKbyAQMs0UGtMyIiIiqAE9sff1it8dlnEB8PzZvDhhy/ciL4yEWPi4yJo/+PG9gf+x29f+7NvlP7uG9PAL3mONzd8GVY+g5ky4bKXrlm0dFw5gxky2YT2o4dg9GjrdXh75GBIiIiogI4kSxdagfbvvsOgoNts619e5saW7jLxcWvSxxnPAtYde4Llv9wkHtOZObzaVDtpvLw7QgoV847T0JSp9On4aOPbEpb/fpW9FasCJs3a8dXRETkElQA34D4ePjhByt8f/kFsmeHHj2gTRvIlevC40KyBhN+IpJzfhs55ZlBlN8fxPkdpsjJrAyf4VDneABO/7Hw9NO2cydyLY4dsyizYcPs79Wq2Zzs81T8ioiIXJIK4ASIiYEvvrBWh99/h9BQGDoUnn8eMmT47+M73F+cVtPf4ojfOHBccOGNRUG8s/AEfv9rbQeSsmdP/iciqVuPHjB8ODRoYCHSldQwIyIici1UAF+H80OzhgyBPXugdGmYONH6fAMCLv058/6cR/81b3DE8yu4dp8nHpzs2fBb+f1/5xz/ZfracAbO3sy+E5GEZA2mY50SNCqX75KPlTRi924LjW7Z0locunSBl16CMmWu/rkiIiLyN11vvwZHjsBbb0HBgtCunfX1zpgB69fDE09cuvhdtmcZNSfUpPak2oSf3EuXqDsJjrXiN9A/gAeGfH3F4rfrN+sJPxGJC4SfiKTrN+uZvjY8SZ+npFBbt8Jzz0HRojBqlM3OBsifX8WviIhIAmgH+Ap27rRzRWPH2iCLhg2hc+crX2kOOxBG9wXdmbFlBrkz5GZo3md4sc9sgnavpMHLDVj40K1UL/XAFWPNBs7eTGRM3EX3RcbEMXD2Zu0CpzVt2sDIkfYu66WXbGpbgQLeXpWIiEiqpgL4En77zQ62ffWVnUlr2dLqjlKlLv85m49s5q2Fb/HV71+RNSgrfW7vwCsfhZFx5ji47Tb4aiqVKla8plizfScir+t+8THr1sGtt9oPX6FCFify+uuW5ysiIiI3TAXwX1wXFi+2wnfWLMiYEV57zVoe8uf/7+OX7VnGwp0LKZGzBDO2zGBC2ASC/YN5o3JnOiz3kLX5EPD3t2kYbdrY36/R+dSIS90vPmzJEjsQOXs2TJliiQ4dOnh7VSIiIj4nzRfA8fHw7bdW+P76q8WX9e4NrVvbPIFLWbZnGTUn1iQqNgoXlwC/AF67+zW6xFcm9yvdLX/1kUes+M13/S0LHeuUoOs36y9qgwgO8NCxTomEPk1JyRYvhp49Yf58+wHs1w/uu8/bqxIREfFZabYAjoqCTz+1Q/WbN0ORIjBihEXxBl9ho/VY5DG6zuvKudhzADg4vF62Ff0mHIAvHrGDSj/+CHXqJHht5/t8lQKRBsTFwbPPWsTI4ME2QeVSWXoiIiKSaNJcARwRAWPG2Obsvn02dO3LL+1q85W6FE5Hn2bo8qEMXDqQiKgIPI4NGQh0/WjYZRzsiLWoiC5dEmXsbKNy+VTw+iLXhXnzLL/3s88gfXq7BFGkyJXfeYmIiEiiSTMF8MGD8P77tst78iTUqgXjx0Pt2uA4l/+8c7HnGL1qNO/+/C6Hzx6mUclGVM/bhp+m/UiRLZ/QfP0xiuavBN9/DMWKJdvzkVTGda239513YNkyayzfutUOSJYu7e3ViYiIpCk+XwBv2waDBlmxGx1tO72dO0OFClf+vNj4WMavG0/PRT3ZG7GX2kVq07tGbw7vyMKx9l2Yvup7jmTIyjs1OzO/bDX6nk5Po2R5RpLqHD8OdevCihUWYTZyJDzzDKRL5+2ViYiIpEk+WwCvWWMH26ZMsdaGp5+2A/VX26SNd+OZ/PtkeizowdZjW6mYvyITGk2gZqEaMG0ah5+pT46IY0y8oz6Dqz7BqXQZIDZeGb1yMdeFTZvgllsga1brDX/+eXjqKQgM9PbqRERE0jSfLIB/+AHq14fMmaFTJ4szu1qE6tLdSxm9ejRL9ixh+/HtlM1dlu+af8eDxR/E2b0bGjSAGTM4lLsIzzfsRljIxYkMyugVwArfWbOsH3zDBtixw374Pv/c2ysTERGRv/hkAVy7Ngwdaru+WbJc/fHDVwzn1R9fJd6Nx8Hh7Wpv82a1N/GLjbOT+W+9ZY3CQ4bwv7Nl2H0q+j9fQxm9aZzrwpw59rPy6682wGL4cMiRw9srExERkX/x8/YCkkJgoO36Xq34XbVvFXU+rUObWW2Id+MBcHHYdugcfr+usEbhjh2tot64Edq14/UHShEc4Lno6yijV/jzT6hXD/bvt5iRzZst3iwgwNsrExERkX/xyR3gq/nj8B+8ueBNpv4xlUyB2cgW14ATfj/iEouDPzd/tBh3Xn+ckBCYNg0aXTjepoxe+duiRfDzz9C9+4X852rVdLhNREQkhXNc103Sb1ChQgV31apVSfo9rtWuE7vouagnE8ImkD4gPe0rtWfWsvIcPOlHlPMHhY5+T/eFa6m77QxTKjfm0VnjIFMmby9bUpqff7ZWhwULbNLfxo3WcC4iIiJe5TjOatd1r5L1lUZ2gA+ePkifn/swavUoHBza3t2WrlW6kjN9Tib8OJP8Jw/Se/aXVN+xmt/y3kyDJ9uwMe/NPKriV/5p2zabkT13LuTJY43mrVppgIWIiEgq47MF8LI9y5i1bRa7T+5mysYpnIs9x7PlnuXNqm8SmiXUHhQby+vrv+e5ueOJd/zoWesFJtzxIPF+HvLpUJucd+6cTffLnNmK4MGD4aWXbIqbiIiIpDo+WQAv2LGAOp/WISY+BoDaRWoz/IHhFM9R/MKD1qyBF17glTVrmF/sbt6o/RL7M+cCdKhN/rJxI/ToYTOzlyyB3LmtAPbzybOjIiIiaYZP/iaf9Nukv4tfj+OHWktcAAARMklEQVShZqGaF4rfM2cs2eGuu6yw+fprIr6cgl+BAjhAvqzB9H24rA61pWV//mkDK8qWtfHFtWtDbKx9TMWviIhIqueTO8DPl3uez9d/Tmx8LIGeQKoXqm4fmD3bLl3v3Gm9m/36QbZsNAIa3ZHfiyuWFGPePBtb7O8P7drZ3Oxcuby9KhEREUlEPlkAVy5QmQVPLWDhzoVUL1SdSumKQsuWNo2rZElYvBiqVPH2MiWlOHzYJrbddRfccw+8/roFSYeEeHtlIiIikgR8OwbNdWHCBGjfHk6dgm7doGtX5bSKOXnSDrS9956lOmzZohYHERGRVOxaY9B88rf99LXhNOswiSWFb4dnnuFogaKwbh28/baKX4GzZ2HAAChcGHr1sgluM2ao+BUREUkjfO43/vS14Xw7eCIT33+esvu30a3Oy1Sp/zbTo64yF1nSjtmzrbe3UiVLA5k82VpjREREJE3wuR7ggbM3czx3Mb4uW5sPKj3KoUw5INZl4OzNSnZIq1wXpkyB48ft8GOjRrBiBdx5p7dXJiIiIl7gczvA+05EcjYwmDfvb23F7z/ulzTop5+s0G3WDMaPt2LYcVT8ioiIpGE+VwCHXGaC2+XuFx+1YYPl9953n6U8jB8PP/9sxa+IiIikaT5XAHesU4LgAM9F92myWxpyPtUkMhJ++80SHrZsscEWHs+VP1dERETSBJ/rAT7f5ztw9mb2nYgkJGswHeuUUP+vr9u3D3r2tB3eUaOsxWH3bggK8vbKREREJIXxuQIYrAhWwZtGRERYpNmQITau+OWXL/T5qvgVERGRS/DJAljSiHnzoHlzOHIEWrSwTN8iRby9KhEREUnhfK4HWHyc61qcGUDx4tbqsGoVfPaZil8RERG5JiqAJfX45ReoXBmaNrVCODQUfvgBypf39spEREQkFVEBLCnfpk02vKJKFTvY1qKFt1ckIiIiqZh6gCVlmzkTGjaE9Omhd29o2xYyZPD2qkRERCQVUwEsKc+ZM7BrF5QqBdWqQfv2dsud29srExERER+Q4BYIx3HyOI6zNjEXI2lcXByMGwfFikHjxvbvjBmhf38VvyIiIpJobqQHeBCg+cKSOBYsgAoV4Nln7XDbJ59ocpuIiIgkiQS1QDiOUxM4AxxI3OVImjR3Ltx/PxQoAJ9/Do8+Cn46nykiIiJJ47qrDMdxAoE3gS5XeEwrx3FWOY6z6vDhwzeyPvFVx47BokX291q1bHzxpk3w2GMqfkVERCRJJaTS6AKMcF33xOUe4LruGNd1K7iuWyFXrlwJX534nuhoGDoUbr7Z8nzPnbOC98UXIVgdNSIiIpL0ElIA1wZedhxnIXC74zgfJ+6SxCe5LkyfDqVLQ7t21u87fz4EBXl7ZSIiIpLGXHcPsOu6Vc//3XGcha7rPp+4SxKftGaNJTvccotNb6tbFxzH26sSERGRNOiGmi1d162eSOsQX3TwIEyebH8vXx5mzIDffoN69VT8ioiIiNfotJEkvuhoGDTI8nyfecYOvAHUrw/+mr0iIiIi3qUCWBKP69oub5ky0LEjVK1qrQ/Zs3t7ZSIiIiJ/03acJJ69e+Hhh6FIEevzrVfP2ysSERER+Q/tAMuNOX4cPv4rCCQ0FObNg/XrVfyKiIhIiqUCWBImLs6GVxQrZhm+mzfb/VWqQECAd9cmIiIicgUqgOX6LVwId9wB//uf9fuuXg0lSnh7VSIiIiLXRD3Acn1On7YJbhkzwpQp1vOrSDMRERFJRbQDLFcXFQVjxljbQ8aMMHs2/PEHNGmi4ldERERSHRXAcmWzZlmbw4svwpw5dl/58hAc7N11iYiIiCSQCmC5tD//hIYN4YEHwM8PfvxRyQ4iIiLiE9QDLP/lutbbu20b9O8PbdtCYKC3VyUiIiKSKFQAizk/xa1GDevz/eQTyJMH8uXz9spEREREEpVaIAS2bLFWhwYNLNsXLOZMxa+IiIj4IBXAadmZM9Ctmx1yW7oUhg6F117z9qpEREREkpRaINKyF1+Ezz6DJ5+0Xt+8eb29IhEREZEkpwI4rdm1yw603XQTvPkmtGoFVat6e1UiIiIiyUYtEGlFdDT06we33AKdOtl9JUqo+BUREZE0RzvAacHChdC6tU1ve/hh6NPH2ysSERER8RrtAPu6jz+2aLNz52DmTJg6FUJDvb0qEREREa/RDrAviouDI0csx7dhQ9i3Dzp21PhiEREREbQD7HtWrYKKFS3TNz4ecuWCHj1U/IqIiIj8RQWwrzh1Cl59Fe66C/butfHFjuPtVYmIiIikOGqB8AUbN0KdOhAebofd3n0XsmTx9qpEREREUiQVwKlZfDz4+UGRIrbz27GjtT+IyP/bu7cYKeszjuPfR+RUtZS1eCKF1NgLtYDnCpUIQg2pYg0hwaiJCtoG0FSNFxI0RKPGXkhCoOChvaiNta1pjNHURaNBVGrlUOWgNqQVLVaR2lpQo6vw78U7lGWZgZlhd96Zeb+fZMLszH+yzz55mPntu/+ZV5KkitwC0Yp27YLFi2HMmOx0xoMGZZ/uYPiVJEk6KANwq3n9dRg3Ltvve8IJ2d5fSZIkVc0A3Cq6urIzuJ15JmzZAo88Ap2dcNxxeVcmSZLUUgzAraJ/f1i9Gq65Jjuj2+WX+ykPkiRJdTAAN7Nt22DmzOxEFhGwfDk89BB0dORdmSRJUssyADejlODhh+Hkk7OtDqtWZbcPGJBvXZIkSW3AANxstmyBKVPgqqvglFPgtddg+vS8q5IkSWobBuBmc/fd2RHfJUtg5crsKLAkSZJ6jQG4GbzxRvbGNoB774WNG2Hu3OwkF5IkSepVJqw8dXXBXXfB6afDjTdmtx19NIwcmW9dkiRJbcxTIedl9WqYNQs2bIAZM2DRorwrkiRJKgQDcB46O+Gii7KTWDzxBFxySd4VSZIkFYZbIBppx47s3wkTYP78bO+v4VeSJKmhDMCN8MknMGcOjBkDO3fCoEFw550wZEjelUmSJBWOAbivPf88jBoF998P06bB4e46kSRJypMBuK98/jnMng2TJmVncHvpJbjvPhg8OO/KJEmSCs0A3FcGDIDNm+Hmm7OzuY0bl3dFkiRJwgDcu3buhJtugvfey05i0dnpUV9JkqQmYwDuLc89l+31XbQInn02u839vpIkSU3HAHyodu7M9vpOngwDB2Z7fa++Ou+qJEmSVIEB+FAtWAAPPOBeX0mSpBbh3+jr8emn8NFHMGIE3HYbTJ9u8JUkSWoRHgGu1apVcNpp2Wf67t4NHR2GX0mSpBZiAK7WF1/AvHkwfjx89VX26Q6H2T5JkqRW4xaIarz7LkydCuvXw7XXwsKFcNRReVclSZKkOtQcgCNiCPBboB/wKTAjpdTV24U1lWOOyS5PPgkXX5x3NZIkSToE9fwN/wpgYUrpQuADYErvltQkNm+Gyy6DHTtg0KDss30Nv5IkSS2v5gCcUlqaUiqd6YFhwIe9W1LOUoKlS7M3ui1fDps25V2RJEmSelHd7+KKiLHA0JTSK2Xu+3FErImINdu3bz+kAhtq61aYMgXmzoXzzoONG2Hs2LyrkiRJUi+qKwBHRAewGJhZ7v6U0oMppbNSSmcNGzbsUOprrBtuyM7ktnQpdHbC8OF5VyRJkqReVs+b4AYAjwHzUkrv9H5JDfbxx9DVlb3JbdGi7PpJJ+VdlSRJkvpIPUeAZwFnAPMjYkVEzOjlmhrnhRdg9GiYNSv7esQIw68kSVKbq+dNcMtSSkNTShNKl9/1RWF9qqsLbr0VJk6EgQPh9tvzrkiSJEkNUrwTYbz9NkyfDuvWwXXXZSe1OPLIvKuSJElSgxQvAA8ZArt3w+OPw6WX5l2NJEmSGqzuj0FrKdu2wS23wJdfQkdHdvTX8CtJklRI7R+An3oKRo2CJUtgzZrstoh8a5IkSVJu2jcAf/YZzJkDU6fC8cfD2rWe1EKSJEltHICvvBKWLcu2Prz6Kpx6at4VSZIkqQm075vgFizITmk8aVLelUiSJKmJtG8AHjMm7wokSZLUhNp3C4QkSZJUhgFYkiRJhWIAliRJUqEYgCVJklQoBmBJkiQVigFYkiRJhWIAliRJUqEYgCVJklQoBmBJkiQVigFYkiRJhWIAliRJUqEYgCVJklQoBmBJkiQVigFYkiRJhWIAliRJUqEYgCVJklQoBmBJkiQVigFYkiRJhRIppb79BhHbgXf69JtU9k3gXzl971Zkv2pjv2pjv2pjv2pjv2pjv2pjv2qTZ79GppSGHWxRnwfgPEXEmpTSWXnX0SrsV23sV23sV23sV23sV23sV23sV21aoV9ugZAkSVKhGIAlSZJUKO0egB/Mu4AWY79qY79qY79qY79qY79qY79qY79q0/T9aus9wJIkSVJP7X4EWJIkSdqHAViSJEmF0vIBOCJ+GRF/iojbDmVNEUTEkIh4OiKeiYjHI2JAmTWHR8S7EbGidBmVR63NoNpeRMQdEbE6In7e6BqbTUTM7tav1yLigTJrnDEgIo6NiBdL1/tHxJMR8XJEzDzAY6pa14569GtEaXaej4gHIyIqPGZ4RGztNmsH/WzQdtGjX1X3oaivlz36dUe3Xr0VEfMqPKaQ81UuS1Q7N800Xy0dgCNiGtAvpTQWODEivlPPmgK5AliYUroQ+ACYUmbNaODRlNKE0mVDQytsLgftRUScCZwHnAN8GBGTG11kM0kpLdvTL+BF4KEyywo/YxExFPgVcETpphuAtSml7wPTI+KoCg+tdl1bKdOvnwCzU0oXAN8CKv0S9T3g7m6ztr3vq81fmX5V1Yeivl727FdKaUG357GNwMMVHlrI+WL/LHEZVcxNs81XSwdgYALw+9L1Z8iCSD1rCiGltDSl9Gzpy2HAh2WWnQtcHBGvln5TO7xxFTadanpxPvCHlL2bdDkwvqEVNqmIGA4cm1JaU+ZuZwx2ATOAHaWvJ7D3eWolUOkD5Ktd12726VdKaX5K6c3SfUdT+YxT5wLXRsS6iLin78tsGj3nq9o+TKCYr5c9+wVARJwNbE0pvVfhcYWcrzJZ4kqqm5sJVa5riFYPwEcAewbz38Cxda4plIgYCwxNKb1S5u7VwOSU0jlAf+CHDS2uuVTTC+ervLnAsgr3FX7GUko7Ukr/7XZTtXNUyHkr0y8AImIGsCml9M8KD32a7EX3bGBsRIzuuyqbR5l+VdsH52tfPwUWH+ChhZyvPfZkCeAftODzV6sH4E+AwaXrR1L+56lmTWFERAfZf+hK+wfXp5TeL11fAxTiT2AVVNML56uHiDgMmAisqLDEGdtftXPkvJVExInALcCNB1i2KqW0M6W0C/gLxZ21avvgfJVExDeAY1JKfzvAssLOV48s0ZLPX60+3GvZewh9DLClzjWFUHrT22PAvJTSOxWW/ToixkREP+BS4PWGFdh8qumF87W/8cCfU+UPGXfG9lftHDlv/H/P5qPAzApH7vZYHhHHR8TXgAvJ9nMWUbV9cL72+hHwx4OsKeR8lckSrfn8lVJq2QvwdbIXz4XAm2QNvesga4bkXXeO/ZoN/IfsyNwKYEGZfn0XWA9sINvcn3vdOfZrn14AHcAveqw5DHgZWAT8Ffh23nXnfQHuAaaVrp/ijB2wVytK/44ENpXmaDXQD7gAuL7H+v3W5f0z5NSvnwHvd3suO79CvyYCb5Xm7fpG15v3pVu/9utDhf+bhX693NOv0vXfAGd0+9r52vtz98wSV/Wcm1aYr5Y/E1zpSMAPgJUppQ/qXSPVKyIGAxcB61JKf8+7HrWmiDiB7OjI8nSAo5rVrpPq4eul6lHt3DTTfLV8AJYkSZJq0ep7gCVJkqSaGIAlSZJUKAZgSZIkFYoBWJIkSYViAJYkSVKh/A9sro+VhO9U9QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(x1, y2, 'o',label=\"data\")\n",
    "ax.plot(x1, y_true2, 'b-', label=\"True\")\n",
    "prstd, iv_l, iv_u = wls_prediction_std(res)\n",
    "ax.plot(x1, res.fittedvalues, 'r-', label=\"OLS\")\n",
    "ax.plot(x1, iv_u, 'r--')\n",
    "ax.plot(x1, iv_l, 'r--')\n",
    "ax.plot(x1, resrlm.fittedvalues, 'g.-', label=\"RLM\")\n",
    "ax.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: linear function with linear truth\n",
    "\n",
    "Fit a new OLS model using only the linear term and the constant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.45993387 0.40392682]\n",
      "[0.39334468 0.03389217]\n"
     ]
    }
   ],
   "source": [
    "X2 = X[:,[0,1]]\n",
    "res2 = sm.OLS(y2, X2).fit()\n",
    "print(res2.params)\n",
    "print(res2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate RLM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4.93558725 0.49914381]\n",
      "[0.10844706 0.00934424]\n"
     ]
    }
   ],
   "source": [
    "resrlm2 = sm.RLM(y2, X2).fit()\n",
    "print(resrlm2.params)\n",
    "print(resrlm2.bse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare OLS estimates to the robust estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VNXWwOHfmUmFAKGXkNAJLXSkQwARVFCaIJcivQmKBQFBAQVpFoQgRVARBCvCp6AoIAgaFJAO0ltCLwkkpM6c749FaCakZ1LW+zw8gZnJzB7v3Fln7732WoZpmiillFIq41kcPQCllFIqp9IgrJRSSjmIBmGllFLKQTQIK6WUUg6iQVgppZRyEA3CSimllINoEFZKKaUcRIOwUkop5SAahJVSSikH0SCslFJKOYhTer9AoUKFzNKlS6f3yyillFKZws6dO6+Yplk4KY9N9yBcunRpduzYkd4vo5RSSmUKhmGcTupjdTlaKaWUchANwkoppZSDaBBWSimlHCTd94TjExMTQ1BQEJGRkY54+XTl5uZGyZIlcXZ2dvRQlFJKZXIOCcJBQUHkyZOH0qVLYxiGI4aQLkzT5OrVqwQFBVGmTBlHD0cppVQm55Dl6MjISAoWLJitAjCAYRgULFgwW87wlVJKpT2H7QlntwAcJ7u+L6WUUmkvSyRmrdoVTONpGykzZg2Np21k1a7gNHvuiRMnsmnTpnjv2717N7t3706z11JKKaXu5ZA94eRYtSuYsSv3ERFjAyA4JIKxK/cB0KGWV7q+dlwArlmzZrq+jlJKqZwp0wfhmesO3wnAcSJibMxcdzjFQfj69es888wz2Gw2TNOkbt26tG3blvDwcMqXL8+nn37K2LFj+f777wFYunQpGzZsICwsjC5dutz3OKWUUiqlMv1y9LmQiGTdnhQLFy6kXbt2/Pbbbzg7O3P+/HlGjBjB+vXrOXXqFBcvXmTq1KmMGTOGMWPGsGHDBoB4H6eUUkqlVKYPwiU83ZN1e1KcPHmSGjVqAFC3bl2cnZ1ZtGgRPXr04Nq1a0RExB/gk/o4pZRSWcS1aw59+UwfhEe18cXd2Xrfbe7OVka18U3xc/r4+HDgwAFA9n0XL15Mly5dWLFiBblz5777Ou7u3Lp1C5AzwAk9TimlVBazZw/06wclSsDOnQ4bRqYPwh1qeTG1kx9enu4YgJenO1M7+aUqKWvQoEF89913+Pv7c+PGDVq3bs3UqVNp2bIlAMHBkn3dunVrVq5cSePGjdmyZUuCj1NKKZWFRERA8+bw1VfQvz8UKuSwoRimaabrC9StW9d8sJXhoUOHqFy5crq+riNl9/enlFJZys2b8Nln8OuvsHo1GAb89hvUqAEFCqT5yxmGsdM0zbpJeWymnwkrpZRSKXLqFLz6Knh7wwsvwOXL8gegRYt0CcDJlemPKCmllFLJtmUL+PvLrLdLFxg5Eho0cPSo/kODsFJKqawvOhq++QbsdujVSwLuhAnQt6/MhDMpXY5WSimVdV25AlOmQOnS0LOn7P0CODvDm29m6gAMSQzChmEUNQxjywO3/WAYhtZzVEop5RgBARJkx48HPz/46SdJvspCEg3ChmHkB5YAue+5rQdw3DTNLNvdICAgAH9/f9zd3fH3979TolIppVQmZbfDmjVw9qz8u3Jl6N0b9u+HdeugbVuwZK0F3qSM1gZ0A24AGIZRAHgPuG4YRot0HFu6Gj58OJs2bcLLy4tNmzbRsWNHRw9JKaVUfMLCYO5cqFQJ2rWDRYvk9latYMECqFrVseNLhUQTs0zTjAu+cTe9BHwDLACmGoaRxzTN/0vpAEaOhLTuFlizJsyalfzf8/f3p169euzdu5d169YxceJE/P398ff357Pb+wxdu3ald+/eXLp0CT8/P+bOnZu2g1dKKXXXuHESgEND4ZFHYMUK6NzZ0aNKMymZt9cC5pqmeQH4GvB/8AGGYQwyDGOHYRg7LsedycoCtm3bRsOGDVm3bl2Cj1m4cCHVqlXj999/5/z58+zduzcDR6iUUtmcaUpJyTjnzskyc2Ag/PUXPPusJF1lEyk5onQMKAv8C9QFTj/4ANM0FwILQSpmPezJUjJjTS/VqlWjU6dO8d4XERGBu7s7hw8f5s8//2TTpk2EhIQQHBxM9erVM3ikSimVzURHw9dfS1DYuVMCcfXq8MknctY3m0rJTHgGMNwwjD+AZsAnaTskx/Hw8Ljv3y4uLsTN5H/++WcAfH19GTlyJJs2bWLy5Mn4+Phk+DiVUirbuHEDJk+GUqXkfG94OMybB+XKyf3ZOABDMmbCpmn63/55DngivQaUmTz11FMMGzaMDRs2ULBgQQAGDhxI3759+fTTT8mbNy/Lly938CiVUioLCgsDDw9Zfp4xAxo3liSh1q2zXIZzamgDh3SQ3d+fUkqliM0GP/4IH34ofXx37ZKZ7pUrDu1klNa0gYNSSqnM48YNCbwVK0KHDnDsGPzvfxAbK/dnowCcXFo7WimlVPowTZnprlolS82NG8O0adCxIzhp+AENwkoppdKSacKmTZLl3KKFBN9u3aS6Vb16jh5dpqPL0UoppVIvMlKOE9WsCS1bwp9/gtUq97m6agBOgM6ElVJKpV6PHrBypTRSWLwYuncHd3dHjyrTy5FBOCwsjF69enH58mXKlSuHt7c3lSpVomfPnnceEx4eTs+ePbl27Ro+Pj58/vnn95buVEqpnG37dpg9W/Z4vbzgtddg+HDw98/2Z3vTUo5cjp4zZw4VKlRg69atREVF8fXXX//nMUuXLqVhw4Zs3rwZV1dXHjxmpZRSOU5MjFS1atRI6jivXn23+H/9+rIHrAE4WRw/E3ZAB4e//vqLAQMGANCkSROKFy/+n8d4eXmxZMkSOnbsyKK4jh1KKZVTRUZKctWpU1LNatYs6NsX8uZ19MiyNMcHYQe4efMmuXNLe+RcuXJx48aN/zymffv2RERE0KlTJ1q0aMEHH3yANS7JQCmlcoIDB2D9enjxRXBzg4EDZc/3iSfuJl2pVHF8EHZAB4e8efMSFhYGyN5v3niu5I4ePUrbtm3p3LkzPXv2ZNmyZTz33HMZPVSllMpYdjusXSvFNdavl+Sq7t2hSBF4/XVHjy7byZF7wvXr12fTpk0AbNmyhbVr1/7nMYsWLeL777/HarVSrVo1IiMjM3iUSimVwXbuBF9faN8eDh2CKVPgzBkJwCpd5Mja0XHZ0RcvXqRChQp4e3uzfPlyChQoAECfPn3o1KkTPXr0wDRN8uXLx4oVK8iVK1eSnt/R708ppZLs+HGp41yvnvzs3BkGD5af2ahvb0ZKTu3oHBmE01t2f39KqSzONOG332TJ+YcfoEEDKa6h0kRygrDj94SVUkplnNWr4Y03YN8+aZwwbhwMHeroUTnEql3BzFx3mHMhEZTwdGdUG1861PLK0DFoEFZKqewuKAg8PaV/76VL0q/3k08k4crNzdGjc4hVu4IZu3IfETE2TBOCQyIYu3IfQIYG4hyZmKWUUtmeacoSc7duULo0fPqp3N6vn/Tx7ds3ywfgVbuCaTxtI2XGrKHxtI2s2hWc5N+due4w4TcNrv9ekUtfPYJpQkSMjZnrDqfjiP9LZ8JKKZWdmCZ88YWUlNy+HfLlk6JI7dvL/dnkfO+9M1lI3kw2PBwO/lyCG9vKYS+8HadGs4i0l8fdWpFzIRHpPvZ76UxYKaWyg/Bw+WkYsHAh3LgBc+fKUvS778psOBuZue7wnQAcJ7GZbFQUBARIwa+QzZWwtJoB/ZsQW/kjLuceQ5TlECU8M7bpRI4Mwn369KFWrVo0bNiQZ555hh49erB169b7HuPv789rr70GQIMGDZg4caIDRqqUUon45x/o0weKF4eLF+W2776Dgwdh2DDZB86GEpqxxnd7bCx89pkcgR4xAoo/EkjZqU2JfWQ8YIIBJrHYnA4wqo1v+g78AVkmCAeeDWTqlqkEng1Mk+ebM2cOgYGBeHh4sH79+ngfs2fPHux2O/v370+T11RKqTQRGwvffgtNm0KdOvL33r2l2hVA4cKSfJWNJTRjvfd2u13+0/j5yRZ47rJ7qP9he3bXaUSY0xHalR2AxXAF04LFcObl5h1zXnb0yJ9HsvvCwxs4hEaFsvfiXuymHYthoXrR6uRzzZfg42sWq8mstomXwzRNk7CwMFxcXOK9PzY2lmPHjuHj45PocymlVLozTVluPnsWunaVJeb335dkq3wJfydmR6Pa+N63Jwzg7mxlVBtfTBN+/hnGj4d/LgVSoOm3VO+/l73h6/GM8GRKyym8UP8FPFw8CDzbj02nNuFf2p+G3g0z/H04PAgnRWhkKHZTrvDspp3QyNCHBuGkGDFiBNeuXaN9+/a0bNky3seUKFGC9evXU6dOnVS9llJKpcr+/TBnDoSGwpdfQpkykvlcr162SbRKrrgZ64PnfAuGedGsGWzdCkWafY/lqWe4ho1r4dC7em9mtZ1Ffvf8d56noXdDhwTfOA4PwkmZsQaeDaTV562ItkXjYnXhi05fpPo/2pw5c9i6dSuurq7s2bMn3sfUqlWLzz77jO7duxMaGpqq11NKqWSx2WDNGsly3rBBjhPFLTlbLFLlKofrUMvrTjDeuRPGvy4z4CJlL9Ji+lR+jwzAbspM2WpYqVSo0n0BODPIEpsGDb0bsqH3Bt5u8TYbem9Is6uWwYMHs3jxYmw2W7z3165dmz179lCtWrU0eT2llEqy2bPh6afh8GGYOlWynBcsyPZ7vcl16BA88wzUrQvb9lyn5ZRxhPUvy++RATxe4XHcnNywGlZcrC74l/Z39HD/w+Ez4aRKjyWD/Pnz07JlSxYtWsT27dvxuJ1F+Prtdl21a9emevXqOGsRc6VUejt8WM7PPPqoBN8ePcDLCzp21EYK8Th1CiZNgiUbA3HyXUedCRc45vYVG6NCeLbas0zyn0TFghUJPBvo0D3fxGgDh3SQ3d+fUiqN2O3wyy/SSOHnn8HFBd56C0aPdvTIUi296jKfPy8dFhcuBLPUZuw9WmM3YgBo7N2YuU/MpUaxGql+ndTQBg5KKZUVdOggXYyKFZPgO2gQFC3q6FGlWmqqWSXk2jWYMUNW6aNjY2ny/BL+KfwKN2MkAFsMC09WeNLhATi5dHNBKaUyyokTMsuNq27Vpw8sWwanT0tno2wQgCFl1awScvMmTJ4sCeHTZ9ip3fsrfKZVZbPnAErmK4Gr1RWrYcXV6pop93wT47CZsGmaGIbhqJdPN+m9vK+UymJMEzZulCncDz/IkaLWrWXvt1MnR48uXSSnmlVCIiNh3jx45x24csWkQe81XK81nj9C91DNvRqrn1xN+4rt2Ra0LVPv+SbGIUHYzc2Nq1evUrBgwWwViE3T5OrVq7hl8c4kSqk0cv06NGkiJSQLF5bqEUOGQIkSjh5Zuirh6U5wPAE3KXWZY2KkxORbb0EQgRTv/AllyweyLfwA5a3lWd5pOd2qdcNiyEKuo8/5ppZDgnDJkiUJCgri8uXLjnj5dOXm5kbJkiUdPQyllKOcOgU7dkCXLpA/PzRsCK+9Ji0Fc8gF+sOqWSXEbpc6JG++CcePQ5F2C6DOMM4bdgiHjuWH8dWzs3C2Zq9McYcEYWdnZ8qUKeOIl1ZKqbRnmrB5syw5r14NuXPDE09ArlywaJGjR5fhEqpmFV9SlmnKKv348bBvH1Roso+KfUdzJPYniNvdMy1sPXyLNXsvZXht5/Sm2dFKKZUaf/wBQ4dKBClYEMaMkX/nyuXokTnUvdWsErJxI7z+Ovz1F5SqdYzG703gz5srsMTmwiO2LeHWjZhmLAZOWGOrMnPdYQ3CSimV4505I81pK1SAQoUk2WrxYujeHdwzth9tVvTXXzBuHGw4HIhHvVXUnPYv+6LWcDnSldGNR/PF+hpYyIOHrRWRln242f1wtVdOVmJXVpGkI0qGYRQ1DGPLA7dVMwzj1/QZllJKZTKmCb//Lnu9ZcrA2LFyu68v7NolnYw0AD/Uvn1yNLpBA9geshbLgGaE1ZjB7sj/o2Oljhx/4ThTH52Kt2cRAFztlckX2xVXuxQ/SkpiV1aTaBA2DCM/sATIfc9tBvA+kL12yJVSKj5ffw21akHz5vDbb/Dqq9JCUCXJsWPQsyfUqAG/BYbgP+kNIp7ugJ1YQJor1C5em2IexQBJ7HJ3vr87VGKJXVlVUmbCNqAbcOOe2/oCv6XLiJRSKjMICpJORiCtBO12+Phj6eU7fTpon/FEBQfLiazKleG7H8Jp9vo0LC+VZZM5maalmibYXKFDLS+mdvLDy9MdA/DydGdqJ79stx8MyagdbRjGJtM0/Q3DKAh8A7QBfjVN0z+exw4CBgH4+PjUOX36dNqNWCml0otpSqLV7NmwciV89500U4iKkrrO2aiuQXq6fBmmTYM5qwKx+azH75HrnCu4nMsRF3mywpNMbjmZmsVqZvrmCimV3rWjpwFjTdOMSajQhmmaC4GFIA0cUvAaSimVcWJipHzk7Nmwezd4esLLL8sSNICrq2PHl0WEhsoq/fvvQ1iBLVj6tsJuxLAHqJmvJqu6r6SRd6M7j8/qhTbSQkqCcHOgwu0AXNMwjMmmaY5P22EppVQGuHVLjhJZLPD22/L3BQukjWDu3In/vgLkP+PcuTL7vXbdToN+3/JvuecJib7bXOGZKs/cF4CVSHYQNk2zYtzfby9RawBWSmUd9y45b90qTRXc3OTvxYvrknMyREfLyay334bz503qdF9LoXrj2HZjD2VyleGW7QY2uw0XqwstSrdw9HAzpSQH4fj2fuO7TSmlMqWICKmLOGeOHCny9ISBA6VTgJtbtq/nnJZsNvjiC5g4EU6ehGrtN1GgxevsvBFIWWtZlnZcSvdq3fk7+O9sueeblrRYh1IqezNNmd3u2CFneatW1SXnFDJN+P576bp48EYghVouo+LQ7ey/tZ0SZgnmPzmffrX63anvrHu+idMgrJTKfkxTlpdnz5YZ7ocfSjejP/6Qhgq65Jwspgm//ipVrnbsgOJtlmE0fI4r2LlyC4bXG86M1jNwd85+xTTSW5IqZimlVJYQEQGffCJZzc2awYYNUKCA3GcY0KhRjgjAq3YF03jaRsqMWUPjaRtZtSs4xc/1xx/QogW0aQPnIo/TcGYvzjfshYkdkEIbJfKU0ACcQhqElVLZx+jR0L//3cIaQUEwYYKjR5WhVu0KZuzKfQSHRGACwSERjF25L9mBePduaNdOFhAOnA2m8bQhXOpaid1R39HDr0eChTZU8iS5WEdK1a1b19yxY0e6voZSKgeKq+U8Zw6MGgX160t9xKAgKS+ZA2a88Wk8bSPB8TQ68PJ0548xLRP9/cOH5brlqz8Cca2+Bt8GxzlsfI/dtDOw9kDGNRtHiTwlskehjVOnYN482LsXfvopzZ42vYt1KKWU49y6BcuXy37vvn2y3NyliwTh8uXlTw6WUKehxDoQnT4Nb70Fn30G1orroF87ooxY9tqhdpGWfPvsIsrkv9sHPssmXZmmbFMEBEgjY8OQqmhxZ8YzmAZhpVTWYbdLF4Bjx6B6dVi0SNoH5vDevfcq4eke70w4oQ5EFy/CO+/A/PlgOt3Cd+h7/FtwMlikuQKmhbPnfdhzyoUy+dNz5Ons5k1YskSqivz7r7SgHDNGilt7eztsWLonrJTKvExTOr8PGyYB2GKRw6mbN8umZf/+GoAfkNQORNevS7Zz2bIQMC+ausPm4jmhHIcKv4mzWQZMZzAtGDhhja3KzHWHM/JtpJ1//4URI8DLS37mzQuffy6NOKZMcWgABp0JK6Uyo7AwqeU8Zw4cPCizlpdflqXmHj0cPbpMLa7T0Mx1hzkXEkEJT3dGtfG9c3tYmKzkT10WSFihjVTpF8mNUsv4M/wUTQs3xXr0JVztVYmyHCLSsg83ux+u9sqJLmdnKjYb/PijLDmvXy/NN559Fp5/Hh55xNGju48GYaVU5rJ3rxwvCg2Vo0affipfoG5ujh5ZltGhltd/2v5FRUmNkilT4JLrHxh9W4IlmoOAr5svP3f4mcfKPUaT6b8RHBKBq70yrvbKd34/oeXsTOXKFamjOW+ebHKXLClveMAAKFLE0aOLly5HK6Ucy26HX36Br7+Wf1epIvu8W7fCzp3Qp48G4FSIjZWj0xUrwosvmhRr+hOFhnbFtEQD0lyhd43etCnfBsMwkrycnan88w/07StBd8wYWWP/7jupqfn665k2AIPOhJVSjhKXKBMQIOdiataEZ54BJyeZyahUsdvh22+lxOSRI1C5zRaqjXydvTe2UtylOC6xLtjM/zZXSGw5O9OIjpY3GBAAgYFSgrRfP1lyrlrV0aNLMg3CSqmMt3gxvPSSBOJHHoGlSyUA59CzvWnJNOXI67hxkrtWrslOag8ezz83f6a4vTjznpxHv1r92HluZ4LnfONbzs40goNlXX3hQkntrlABZs2C556TphxZjAZhpVT6s9slMlSrBqVKyXLh009LtmomS5TJyjZvltXXP/+Ewi2X4zXlHY7HHKBgbEFmtp7J8/Wev1NeMkud842rBR4QACtXSuLVk0/C8OHQurVkzWdRGoSVUuknJEQSq+bOhePHZXo2ebIUI26h/WXTyo4d8p/2l1+gaKWTVH3neQ5E/wQx4GRx4ssuX/Jo2UcdPczku3VLeiYGBEjCXv78MHIkDB0qF3LZQNa9fFBKZW4vvSSJMi+/DMWKSS/fHFbHOb0dPAidO0O9erD933M0mjKMa//z5XDsrxjI0r5pmmwP3u7gkSbT8ePw6qtytnfQILktrhb4zJnZJgCDzoSVUmnFZpNaznEz3IgI6NpVlpxr1XLs2LKQVbuCE02KOnECJk2SrfRcha/QaOJ0/rEG8LctloG1B9K2fFue/fZZom3RWafBQlyWfEAArF0LVit06iRLzk2aZNt8AW3goJRKnbizmR99BGfOwK5dkulsmtn2izO9xHVAioix3bnN3dnK1E5+dKjlxblzspq/cG0glP+ZKg2COZnra8JjwuhVoxcTmk+gbH6ZJWaZBguhoXe3LI4dg6JFYfBgmQF7ZdLksEQkp4GDBmGlVMpcuCBZQCtWQGSkzIBHjID27eWYkUq2hDogFXHOg39UM+bMgegSm6DXY9iNGAD8S/kz98m5VClcJYNHm0r790vgXboUwsOhYUP5/HTuLBWusjDtoqSUSh8xMVJzt2xZOZf5889STGP48Cx1NjOzerA0pD3Kyo0dZTjzd1m226KpP2QxB4qPJixWArDVsPJYuceyTgCOjYXVq2XJedMmcHWVwizDh0OdOo4enUNoEFZKJe7CBTmXOX++tA7ctw/y5JHSgM7Ojh5dthHXAckeYyFsVylCt5XDHulE7rbzyd/yPf4KP0n1gtU5fOUwsfbYrLPfe+mSdLyaN0+Sq0qVgunTpbhGoUKOHp1DaRBWSiVs926YMUMqE8XEQNu2smQYRwNwmnqppS/DJ4Zw+dhl7EVXYW0VjnP15YS7nKJinlosfHotbcu3ZVvQtqyx3/v33zLr/eorqXDVurUsQT/5pCReKQ3CSqkHREZKpnPu3HDgAKxZI60En39eqhOpNGezxZ3g8uJi1Gno2wYs0dgMyO/uxbx239CpcicshpwqzdSFNiIjpQ54QABs3y4rJoMHy2eoUiVHjy7T0XPCSilx5gyMHStnewMC5LZnnpEygbNmaQBOB6YpW6Q1a0LPnkCprRQd1gOs0WCABQsjGwylS5UudwJwpnX2rCTqeXtLCcmbN+VzFBwsvRM1AMdLZ8JK5XS//SZ9e1evln8//bScywTJUs3imaqZ1YYNErP+/ht86u+i1szx7ApfSwFrAZwtzthNOy5WF1qWaenooSbMNCXBKiAAVq2S2556ShKtWrbUI2pJoEFYqZwoOvpucJ0+Xeoejhol5QBLlXLs2LK56UsuMfVtJ0KPF8C59B7KTZjAcWM1N235mf7odIY/Mpw9F/Zk7j3fsDA5WhQQIGW7ChbUz08K6TlhpXKSI0ekqMbSpVJUw8dHlhELFQL3LNC0PQvbuxcGvBDB9s3uUONLrI++jc3jEBbc6OI7mIUdJ5LPLZ+jh/lwR45IYtVnn8GNG1C7tiTqdeumn5976DlhpdRdNpuc5w0IkJ/OzrLXGxsr93t7O3Z82dzRo1Iy+8svwShwBaehA4gt8gtSE8tCwejRBJ1umnkDsM0mHbACAmDdOvn8dO0qS8716+uScyppEFYqu4orG3npEnToAIULS8HhQYOkoYJKV0FB8NZb8Mkn4JLvGg3Gz2CbZRaxRMsDDMCEGMsJzoUkadKUsa5dk8F/9BGcPAklSsgbGjRISkuqNKFBWKnsZs8embVcuAA//ADFi0uj2Xr1cvy53qQ0R0ity5dh6lSJXXbnm9R7ZRYHPN9lW/RNCllaYUbU5arL+5hmLAZOuNn9KOGZiZZy9+yRRL0vvpDjRs2aSd5Ahw45/vOTHjQIK5UdxMTA999L8N2yRfbnevaUJWcnJ2jUyNEjdLgHmyMEh0QwduU+gDQJxCEh8N57cporrNAmSo6cSWjeP9gWE0rHsh15q8VbHAvOz9iV+3CKLkikZR9udj88rdUY1cY31a+fKjExsHKlfH62bpXPT69ecja8Rg3Hji2b0yCsVHawYIEkyJQtK5Ggb19pgK7umLnu8H3diQAiYmzMXHc4yUE4vpl064peBATIZPF6aAy+Q8dzuPBMgjCxxFpY1H4R/Wv3B6BakbixuHAupHK6zcaT7N5ypOfP6+fHATQIK5XVmCYEBsqs5YknZMbboweUKQOPPw6WTF7UwUEebI6Q2O0PenAmHXQliiFjQ4neWYzrV6Fmzy/J7TeBw7eO3/kdA4NL4Zfue54OtbwcF3Th/s9PXDnSxx+X2s5t2+rnJ4Ppf22lsoqICEmUqVMHGjeWxufXrsl9+fNLPV79Ak1QQvuuSd2PjZtJm3aDsH0lCf64ORfXVSayygrKzqjJ7nI9KZQ3DzNbz8TdyR2rYc1cDRYiIqRvb926dz8/w4bJsaO1a+WCTj8/GS5JM2HDMIoC35qm2dQwDB/gc8AOHAMGm+l92FgpBe3awcaNUK18vnYeAAAgAElEQVSaLB/26AEeHo4eVZYxqo3vfTNZAHdn6337sQ9L3Aq+HkH44WJcPxqCrcAarI9YcKqzlAj3gzi5VeSrx7+6U16ysXfjzFNs49Qp6V60aJFctOnnJ1NJtFiHYRj5gRVAEdM0axuGMQVYZprmIcMwfgJGm6a5N6Hf12IdSqWA3Q6//ir7dYsXg6enBGCrVbJV9WxmijwsyD643AwSpN/p6If7JS+69L9BuOUA9GkB1iip7Wz3pIzLQP4d+w5Olky0u2eaUhdzzhzJkLdYoGNHOdurn590l9bFOmxAN2A1gGma4+65ryBwJdkjVErFLyQEliyRqkRHj0KRInDoEDRsKLV4Vao8bD82vsSt6yfy0qdzbm6cggJ+R4h88n/YnKLkTtOggNmed594MfME4Bs34PPPZb/38GH5/IwbJ12MSpZ09OhUPBL95JimeQPAeODKyTCMbsAB0zTPpc/QlMphLl6EcuUgPFyC7sSJ0LkzuLo6emQ5wr0JWlEX8hKyxZfIE0WweO+l9tTB/BP1NW7W3NhtVkzTxGI480rzLv8J6hlxFvk/Dh2SC7clS6Suc/36Eoy7dtXPTyaXoss3wzDKAq8CjyZw/yBgEICPj0+KB6dUthYTI51njh6VdjpFi8Ibb0jj89q1HT26HKeEpzunjlsI2eLLrcPFMYocxWXwEKKLr+aw3Z3xTcfzSqNXOHT5UIL7vel9Fvk+sbHw448y692wQRpydO8uZ3vr1Uvb11LpJskNHAzD2GSapv/tPeKfgQGmae5L7Pd0T1ipB8SdzVywAM6dg4oVYd8+bRnoQKdOQf8Xwtn4Yy6o9H9YH5uIzXMfBgbtyvVhUacpFMldJNHnaTxtI8HxHHny8nTnjzFptJ1w5YokWc2bJz2gvb2le9GAAVKaVDlcejdwGAP4AHNuL1FPME1zcwqeR6mcZ8UKaXgeEwNt2kggfvxxSbhSGe7CBZgyRf5nMHJFUfKlAQTl+VKaKxgWRtdbwNQn+if5+VJ7Fvmhdu6URKsvv4SoKGjRQspztW8vVdFUlpTk/+VM0/S//XM0MDq9BqRUtnLrlgReX19o0kT2eocNkz8VKzp6dDnWtWswcyZ8+CFEcZM6L3zIvwVmEhRz485jrIZB3jyXHvIs/1XC0z3emXCKa0NHRUlBjYAA2LYNcueGfv1kyblq1ZQ9p8pU9GS2Uunh+HF45RXJSB0wQAIxQOnSMnvRAOwQYWEy8y1bFqa9G0mVfrPwfLMc2/O8QYuy/nze4fNUFdoY1cYXd+f7VzUePIucJMHBkh/g4yMV0a5elSuG4GDpDKEBONvQNQyl0tqgQbJnZ7VCp05yNrNJE0ePKkeLjJT6FO+8A5fdt1CszwdYim5lZ/RlWpVoxeSWk2lQsgEA5QuUT3GhjbjkqxRlR5umNN8ICJBmCna7VEEbPlyS9bSaVbaU5MSslNLELJXtXbsGy5bBkCGSXDV/vmw2DhokPViVw8TEyKmdSZMgKNhOmefe5lTpSZiYWAwLs9rOYsQjIxw7yPBwaRsYECAJevnzQ//+kmxVtqxjx6ZSJL0Ts5RSAP/8I2czly+XqZavryRbDRni6JHleHY7fP01vPkmHD1qUrH9D5R5YTwnb9090GFgEBYV5rhBHj8uS8uffCJFWmrUkBWU7t0hVy7HjUtlKF3fUCq5rlyR/rx16kimau/esHu3BGDlUKYpVRpr1ZJYFuO9gQozGnCkztM4uUUyyX+SY5sr2O3w00+yzFyhAsyeLZ2Ltm6FXbtkBqwBOEfRmbBSSXH2LOzdK1+eBQvKecwPPoA+faSus0p3iVWi+u03GDEjkANhmyhcohBVe3zJgYiNeDt5s6j9Ip6r+RxOFidal22d8c0VQkKkg9FHH8GxY1CsmEzTBw+G4sUzZgwqU9I9YaUSYprSNGHuXFi9GvLlk71eLaqR4RJqrjC1kx8lYrwYNw7W/xsozRVu13bO75afCc0nMLjuYNyc3Bwz8H37ZK932TI5rta4sSRadeqkn6NsTPeElUqtDRvky/Lff2XmO2qUJMroF6dDxNdcIfScO/17unDtIOQvf4Tig4Zw3i4B2MDghfov8GKDFzN+sDExctEWEACbN4Obm7QNfP55WSdX6h4ahJWKs38/uLtLE4V8+SBvXkmt7dpVvkhVqqSmscG9FadirucidGtFwg+WgIInqTPpdXbzKbcMJ5wsTpimiYvVhTblMniP/uJF+PhjyY4PDpYz4TNmSHGNggUzdiwqy9AgrHK2uCYKc+fKrKVfP+nfW7cu/PWXo0eXbaS2sUEJT3dOn7UT+mcFwvZ6g8dFXPr0JabUcvZZDIbXHc7YJmM5cf1Exu73mib8/beUk/z6a/k8PfaY7P0++aSWI1WJ0iCscq7334f33pMmCqVLw/TpEoRVmotvOTkixsbMdYfvBOGEZspXrkDhg/UIXJELs9R6rEOfxVboH6KJ5dFSXVncaTo++aRbW1GPohkTfCMj4auvZMl5xw7Ik0e2K4YNk6NqSiWRBmGVc5im1N9t0AAMA4KC5GymNlFId4k1NohvpvzaioN8OT8Pa5bnJTzGQvHBgzlXaDE2A8DCy3UCeK/d0Ax6B7edOSPLzR9/LEfVKleWVZRevSQQK5VMGoRV9nfzpmSnfvSR7Ptu2AAtW8K772opwAySWGODe2fK9hgLN/8pzY1t5TgaY1Kj/2yCykzhXNTdZgpWw6BQvpD7nis1e84PZZpy/ikgQBKuAJ56CkaMkE5G0k1OqRTRbyCVfYWEyBell5csEzo7S0WiBlIjWANwxkmsscG5kAhMm8HNXT6cW9iCkN8rYGn6EdaxpdhT4kWqF6/KgnYLEiy0ETeTDg6JwOTunvOqXcEpH/TNm3ebJbRqBb//Dq+9BidOwPffy4WcBmCVSjoTVtlLbKx8SVasKJWHfvhBZi3Dh0P9+vql6SAPa2xgs4HTiTKc+qU0sXl3YenYA4vPH8Q6X8TDqMSqnl/QqmwrAPyK+MWbeJWUPeckO3JElpg/+wxu3JDKaJ99Bt26aZa8SnMahFX2cP687NMtXCgz3BMn5EzvkSN6tjeT6FDL676AaJoyoRw/Ho4drIz1scnQcAJ2wwTToJCtHws7TKJV2ZJ3fqehd8N4E68S23NOlM0Ga9fKkvMvv8iqSdeuevGm0p2ux6msbe9eePZZ6bs6YYIsHQYE3P3S1ACc6ZimxLlHHpHCUTcL/UbF6Y2xNXoTjNsV/AyDttWK0LF2yYc/2W1xe8tJvf2Oa9ckN6BCBVkxOXAA3n5bypQuW3Y3iU+pdKJBWGU9N2/C9evy96AgWLdO9n4PH5a/P/WUZjpnUn/+KblMbdpAkPk3Vaa15mzLloQ7nWF049F39nzdnVwZ1vCpJD9vYnvO/7F7NwwYIPkCo0bJRdw338DJkzI1L1o0NW9TqSTT5WiVdRw8CPPmSRWrYcNg2jTpQBMUBLlzO3p06iF275bYtmZPIO6PfEnZt3dxwrYFm6UwH7T5gCF1h+Dm5MbTvk+nqNjGw/ac74iOhpUrZaXkjz8kZ+C556ScpJ9fWr9lpZJEGziozG/1avjwQzkm4uIie3UjRsh6pkoT6XW858gRaRb01VeQq+43RLTrjokkUA2qPYj32ryHh4tHql/noc6dk1yBBQukAUe5chJ4+/SB/PnT97VVjqQNHFTWd/Xq3Xq7X34pDdCnTpWKVkWKOHZs2UxqS0rG58wZeOstSSp2KRRErTffZrflY0zkot9qWCntWTr9ArBpytp3QAB8+61kzT/xhCRatWmjx9NUpqGfRJV5mKbUb+7WTfqt7t8vt8+dK9nOY8ZoAE4HDzvek1wXL8LIkZLn9Pl3l6k+6mXsw8uz3+lTOlfujJuTW7znfNPMrVtS+7t2bWjSBH76SVZNjh6FNWukMpoGYJWJ6ExYOd6tWzJl+ugjyU719IQXXpCfAAUKOHR42V2qj/cAX/x+jjETowneWhLT6QZl+k/nktdH7LHd4jm/53iz+ZuU9ixN4NnA9GmwcOKE5AssXixJe35+svzco4fmC6hMTYOwcpywMPDwkM4zr70mhe8XL5YjR7lyOXp0OUZiJSUfJjwchowJZfmiwtiLbcXSty9msb85ablF4yLtWdRxBpUKVbrz+ITO+aaI3Q7r18uS848/ygy3UyfZ723WTI8WqSxBg7DKWNHRd1sH3rwJO3dK796DB8HbW784HWBUG9/79oQhkeM9QFSU5DpNmQIXr7hh7TwCKi/AbgCmQcGokZhX298XgNNMaKhkyM+dK5lfRYrAuHEweDCUTNq5YqXuCA936GqJbo6ojBEcLGmypUrJnu/Zs9C9u1QqAjmnqQHYITrU8mJqJz+8PN0xAC9Pd6Z28os3KSs2Fj75RKqCvjAyFk//T7GMKY2tyoJ7HmVgs1xL1nJ2khw8KLNcLy948UXZpli2TLLA3n5bA7BKnn/+gYED5SJu3z6HDUNnwir9mKYsNbu4wK+/wuTJkqH6/POaoZrJPFhS8kF2uyQZv/EGHDlqp1y77yj5/BscjjhMbsMXl6h2hDh/jmnGYuCEm90vScvZiYqNlfrfAQGwcSO4usp2xfDhUDdJJ0CUuss05WI/PByaN5dJwP/+59DtLw3CKu2FhMhy4bx5skT40ksy+23eHMqUcfToVDKYpiQYjxsHu6/8Sb5WCyne50+ORx+lqkdVZrdfCRH1eP37/bhG+xJp2Yeb3Q9Pa7WHLmcn6vJl6Xg1b56smnh7yxG1/v2hcOG0e4MqZzhwQBL1du+WExi5c8P//R/UqnU3AdRBNAirtLNrl2Q4L18uGc8NGkD58nKfu7sG4Czm99/h9deluFSeR+fA0y8SapiERkPXiiNZ3u1drBYpFWkYBjPXuXAupHLqin3s2CGz3i+/lI3nVq1g9mxo1w6c9OtKJUNUFHz3HcyfD1u2yIrcM8/Id1Pu3FI/NRPQT7VKHZvtbp3mUaOkQEKPHjB0qJzVVOkuratd7dghM99ffoFC1XdQ8s3XCLL8BnHF9UwLvx0K4Yc9F+68TmLL2Q8VFSV1mwMC4K+/5Auyf3/ZtqhSJcXvQ+VQcUvO334LPXtKhbQZM6BvXyhUyNGj+w8NwiplTp6UK8zPP5cM5xIlZOmwUCEtBZiB0rLa1cGDsue7ciXkq3CAam+/wX7b9ziRlzyx7Qiz/nJnz9caWzVlvXrvFRQkS4QLF8KlS5LtNXs29O4tGfNKJVVc7sD8+fDYY/DKK9C5syRdtWqVqfNPNAirpLPZ4OefZcn5p5/kg/3UU7K8A1ImSWWotGhmf/IkTJwISzcF4uT3HeUm7ucEv3Da6sGkppNY/FNlDHKR29b8zp6vq71yyrKfTVPWuQMCpJmw3Q7t20uiVSb/slSZUFCQ9BFftEhqhJcsKWfFAdzcoHVrx44vCTQIq8TFLe+cOydBt0gRmTINHKjHQhwsNdWuzp+XhPWPPwZ8/w/6diLGsHEc+J/f/5jddjYFcxVkXeBGgkMicLVXxtVe+c7vJyv7OTwcvvhCgu++fXK86JVXYMgQzRVQyRP3fQTSjvKXX6Qc6fz58jOL5Q7oZaeKX1wB/F69oGNHuc3bWzILz5yBSZM0AGcCKWlmf/UqjB4tW2ULll6hyshXMZ/pjGnIjNpqWKlWuBoFc0kDjWT36r3XsWPw8stytnfwYMkfWLxYZjDTp2sAVkl36ZJ8Znx95fMDMHOmNHdZs0ZWVLJYAAYNwupBYWGyT1erFjRuLGn8pUrJsiFIUXxnZ8eOUd2RnAB586bUtChbFmbMDqXCwAm4jS7DPo8PaF2udYLNFZJTzAOQz8ratXImvEIFmDNH/v7HH1IgoV8/yZZXKjFx2xfdu8tF/5gxckF3/brc7+eX5S/kktRP2DCMosC3pmk2NQzDGVgJFAAWm6b5ycN+V/sJZxFxSzzvvitZzjVqSIZzjx5S31llWollR0dGwug5gSzasIlbhxtQ9bHtBJWeTmj0NbpU6cJb/m9RuXDl1DdXuH4dPv30bterYsVkuXnQIChePA3fscr24r6PLlyQ4OvhIf2fBw+GypUT/fU46dUnOzHJ6SecaBA2DCM/sAIoYppmbcMwXgbymqY50TCMtUA30zRvJvT7GoQzsbg6zh99JPu7PXrIWuWRI3LGV8tIZmkxMRITx80P5MoTrcApEgz5//vj5R9ncsvJ1C6eBsfI9u6Vvd5lyyAiQlZLhg+XbQwXl9Q/v8o5duyQvd0rV+S7CaRJR6NGya5q9eDJAZBVooeu4qSR5AThpCxH24BuwI3b//YHvr79998BrR2X1Zw9K4lVcXWcT5++m5VasCA0bKgBOAuz26VeSuXKMHiIDXv998E54k4AHlR7EGt7rE1dAI6JkbO9zZrJqsmyZXIRt2uXFEbo1k0DsEqa8HDJE6hbF+rVgxUrJPkzrq78o4+mqKxkWvbJTk+J7mKbpnkDpCLObbmB4Nt/vwYUffB3DMMYBAwC8PHxSYtxqrT01FOwZ4/s0w0bJnWcrdbEf09laqYpRyXHj4d9++2UarsS74FvcjbyEAYGhmHganWlT80+KX+RCxcknXr+fMmWL1NGkmP69dO+zyp54pacFy6U5L0qVSR/oFevNDknnhZ9sjNCSlLJwgB3IBTwuP3v+5imuRBYCLIcnZoBqlS6ehU++0yOh2zeDHnySFGNokWzfEKDumvjRikx+ddfJiWa/0yZd8ZzMuofKntU5tv231LcozibT29O2X6vaUolq4AA+PprmQU/9pgk8D3+uF7AqaSLjpZqMPPny4Vb797w3HMyC27SJE1X4FLTJzsjpSQI7wSaAN8CNYBtaToilXqmCX//LcE2rgZv48Yyi8mTR/Z7VbawbZuUmNx4JJDcjT+jZPttBMXupbR7aZY8voQefj3u1Hdu5NMoeU8eEQFffSXBd+dO+ewMGSLlJH1T0ZxB5TwnT8qM95NP5KhRmTJ3jxMVKABNm6b5S6akT7YjpCQILwHWGobRFKgC/JW2Q1Kp9u+/Emg9PKRe6tChUL26o0el0tDevbKt/3//Bx5NP8UYMIBw7ITHwisNX+GdVu/gYk3hnuzp03IBt2iRrKRUqSLJez17SiBWKinuLarRpYt0MGrXTr6PHnss3aujxSVfOSI7OjmSHIRN0/S//fO0YRitkdnwm6Zp2h76iyr9/fuvfGna7bKnUrmyzGDatoW8eR09OpWGjh6FCRNkgSN36YNUmfQGB82Vd+63GlYKuhe8LwAn6ZiGacqadkCARHaAp5+GESPA318T9VTSXbggiVbLl0NgoHwHzZ8vR9a8vTN0KKlqLJJBUlRexDTNc9zNkFaOEBMjKfzz5sFvv0kBjZ497159du3q6BGqJEpKkDx7VgptfPIJOBc5SdXxEznotIyzzrnpV6UfK/avINoW/Z9CG4k2eLh5U5pwzJ0Lhw5JA47Ro2XZWZMqVVKZpuSczJsne76xsVIL/NIlCcL16jl6hJlWkop1pIaeE04nb7whhX9LlZIvzH79JK1fZSmJnWW8fBn6vxTGD1t3Q/kfcamwB1vxDThbrQyvN5zRTUZTKFehBAttNJ62Md7klIYxV1gRs1OS9m7elMSY4cPlaJGbW0a8dZWdHDok2xb580tRjSFDpCtWDpWmxTpSS4NwGrDbpUj5vHnwwgtyhXnqFBw4IEvOmp2aZSUUJIu6edDG3pyZ79mJLLEO/tcebtd29jAb8WHb2fRrWCfR5y8zZs2dNsAWu40WJ3bw3M4faXZql5zj7dZNgu8jj6Tl21LZ3Y4d8n1ksdzuAIKcj3v0US1JSvKCcNardp2TXLkiJY8WLJAi5UWKyJcmQOnS8kdlaQ+eWbTHWLi5szRn/yrH3+YN3J9+Gyp9KAHYAEwLTrbyLN4cSr8knDQq4elO+PmLdN37K712rcU79CLnPQqy4LF+DF46VVdPVNLduiW5Jh99JEE4Vy5ZgYvbAmvf3tEjzJI0CGdWpgn160sN3mbNZOm5UyetQpTNxJ1lNG0GYXt8CP2zPLYoE7cnpuJafw6hMVdxtfkRZfwLpg0DJ9zsfkkrOLBrF8u3LaTomu9xi41mm3c1prTox9YqjZncpaYGYJU8U6bAO++keVGNnE6DcGYRFiYFNb7/XpZ1nJ3lg+7jA9WqOXp0Kp28/KgvI966ysWjV7EXWYW1zQWsVb4i0ukKzX3acP7M04RG+BBlOUSkZR9udj9c7ZUTLjgQHQ3ffSdZzn/+SalcuTjZoRtveDXnD7dilPB0Z3ImPKahMpmYGFi9WpacX3tNquoNGSI/mzbVbPk0pEHY0fbvlw/60qWSIFOjBgQHy1LzE084enQqnZimJJG+8YYX58NOQp82YInBZkCpPJX5vPN3NCvV7E7iFjGVcbVL95h4Cw6cOyfbFgsWwMWLUL48fPAB9OlDGU9PljngPaosKChIimosWgTnz8skIOx2UURv7ww/YpQTaBB2pG3bpFmCq6scKRo6VLsXZTMPHj969TFfcl32Ytw42LnTxOvR78nbfCg3bDEAWAwLg+r1pFmpZkAiBQdMU3r0BgTI7Ndmkwu34cMzpBiCymZMU2a5p09LOdKFC7UsaQbQ7OiMdPKkzFTy5pVCv3a7/LtrV+lepLKVB48fRQbl58aWSkScyU/RRr+Sq904TkbvoFS+UpwPO4/NbsPF6sKG3hseXt/51i0phBAQII04PD0lQWboUJkBK5UUV69K4ufq1VKoxdlZfpYpo3XlU0mzozMTmw3WrJEl53XrZJbbp4/cZ7HIF6fKluJaqUVfzMv1332JPFEES8UNuL0yhot5dlDKvRSfPv4pPav3ZHvw9njP+d7nxAnJTP3kE7h+Hfz8ZLbyv/9B7twZ++ZU1hTXjGPePMl0joqS2e+FC7LU3LKlo0eY42gQTm8vvSQJVsWLS4GNgQOhZElHj0plgFPHLYRsqcWt8NNQcyLWpwOx5dtNjJmfgMcDGFB7AK5OrgA09G4Yf/C12+HXX2XWu2aNXLh17ixLzmncdUblAIGB0szFwwP695dkKz+/ZD9NkkqhqiTRIJyWTBM2bZKrzHHjJMlqwACpvdu+vSz3qGzv9GmYNAnOfdYcanwBzz4Hhh2bCR6xT1A11/M8/0giSXehoVLNau5cKRhdpIh8poYMAS/9slNJdPCgfB8VKCAfyoYNYckS6Ngxxc04Ei2FqpJFMzfSwvXrMGuWNE5o2RLWr4djx+S+6tXlfK8G4GzvwgUpaFahAnzx42m8X+oNT/cG7LcfYcHdKMKYtjUSfpIDB2SLwssLRo6UXIEvvoAzZ6R4tAZglZjoaFlq9veHqlVly+LqVbnPMKSHbyq6YcVts9wrIsbGzHWHUzHonEtnwqkVEwOVKkmh8rirzGee0dJtOcj16zBzJnz4IUQ6XaDyi1M4kncBlwwLTYo+xR/nfsI0Y7EYzrzcvON/ZwuxsdK5KCBAmnG4ukL37tK3t26ScjuUuuuVV+SzVLo0TJsmSXuFC6fZ0ydUKCZJBWTUf2gQTq7wcFixQrIIv/hCZrizZkkVmRoPmeGobCcsDGbPhneWBhJefC2lh5ziQv7vOGzG0L9Wf8Y3G0/JvCUTbK7A5ctyHnPePGmT5OMjX5r9+0s3I6USE5czMG8evPkm1K4tKymPPy6FNdLheFFclbf4bk8O3VcWGoST6sAB6Yn5+edw44YkM1y5IleY3bs7enQqA0VGysmyd96BS7nXY/R6HCyxnALalG7D3CfmUq5AuTuP/0/S1fbtMlP58ktZOmzVSqbR7duDk/5fUiXB1auSMzBvntSVL1xYkhFq15YJQZUq6fbSo9r4xtv56z8FZB5C95Xv0v/HJ8VPP0kRBBcXOdM7ZAg0aqSZqTlMbKzsNkyaBGfPR1K++zxuVBhPpD0WAKthpXmp5vcF4DuiouDrryX4/v23ZKcOHChLzpUrZ/A7UVlabKzs9V68KBnyb78teSeurhny8g8tIJNED9tX1iCs5DzmggVQtiwMHgwtWsD770vBcl0mzHHsdvjmG1ntO3IshjKdPqVQ7bc4Fh1MvWL12HtxL7H2WFysLviX9r//l8+elc/SwoWy/OzrK0fWeveWoi1KJebWLdkC+/VX+enkJPsglSun6HhRWuhQyytVwVL3le/SIBwnNhZ+/FGWnNetk72U4cPlPjc3Oe+rspXE9qRMU47mvvR+IMdiN1K4RgzF+izjZPRxGhVtxDctl+Ff2v+/e76mCZs3y6x31Sr5d7t2MGKELD3rCopKisOHZbn5s8/kyFq1apIAWrSorMhlYWm1r5wdaBCO06ePJFp5ecGECbJUqMdBsq3E9qQ2bZLKooFn/4S+LcAazWWgfO7y/Nj5R56o8ATG7WB6Z883LEwu4gICJIegQAF49VXZvtDezyo51q+H1q0l8bNLF0m2ykbFWdJiXzm7yJlB2G6XD/n8+ZLZ7OMDw4bJh71dO02OyQES2pN6c9E55h7xYv16k4L1NpC733OEW6JvP8KgftFOPFnxyfuf7OhRKSf56acyY6lVS0pLPvusHlVTSXPuHHz8MRQrJltgTZvCjBmybVG0qKNHl+bSYl85u8hZ0ebyZVnaWbBAMgoLFYJ//5Ug3KiRo0enMtCDe0/Rlz0I2eLL6aPFOFPtT8q9PY7jtk1YzfxgOgF2DJzYur8YqyoG06FGcfjpJy6+8y5F/9xEtMWJTX7NcR35As2fa59tZiwqHcVV2PvoI+kjbrNJhb3BgyXJatQoR48wXaV2Xzm7yDlB+OZN6QwSHi5XmRmcUagyl7g9qZjruQj9owLhB7yg5E7cXuhBaIGNuLoVoUzUcGw3WxFtOUakZR9udj/y3irB2fFvw7+/wokTGB4FeL9JD1bUaMtlj/y4H7cydfc5/XJRiRsyRBL2ChSQnJPBg7ULVg6UfYNwaCgsXQr798uyc548chazQQNJ7xljm3IAACAASURBVFc5Wr9alXn19RhCrpyFKu9h9f8bW+E/sTrnY2qzqYx4ZATV3tyEAbjaK1Pzgiu9/llDxwObcI+VzjNv1O/BCq86xFrv/t8opx6zUEmwe7ckWr3+OpQqJactGjWSJCvdtsixsl8Q3rlTPugrVkhqf716EBEhH/L+/R09OuVgV65IUaq5c4sTVW4l9OkKhg2bAS1KPsPKHgvxdPMEwDuPM35/b6T3Pz9SP+gAEU6urKrSnJ/9u7Dkw4EsG7OG+Lpx58RjFioBUVHw7bey5Pznn3LS4rHHJAg3aSJ/VI6WvYLwkiWS5Zwrl/RYHTxYa+8qQIqcvf++/AnjAr5D3+Go50fYTEnOshpWWlesJQH4wgVYuJB1AR/hfvkiZ/IVZYp/P76u3provJ5M7SRnM/WYhXqoyEgoV06SripUgA8+gOeeg/z5HT0ylYlk7SC8f78kWTVpAt26Sdm/OXNkmSdfPkePTqWzpNSejYiQboDTpsHV8Ov49pvJmeIfctQexZMVnuSXE78QY4uRQhthBaFHD6nMERODe9u2BD7endfCihN0I/o/r6HHLNR97HapMfDHHzB5ssx6X3lFOqm1bCm9oJV6gGGa8S2opZ26deuaO3bsSLsnjIyU5Z358+XD7uoK48fLH5VjPHjOFyQATu3kR4daXkRHwwuTrrP4u73EFvsNF6+TWCqtItJ+g+7VujPJfxIVClYg8NgmNv0YgP+P+2i44YhUserXT85lVqyYpHHoMYsc7upVOZI2f75U2ytaFA4d0hlvDmYYxk7TNJO0DJv1gnCrVtLBqEIFWW7u00d6rqocpfG0jfEvBed1Z5hPS14dG8tF61bo/RhYYsAAd3tV3mnxASP9W8OpU5I7sGgRXLsmyXojRshM2MMj49+QyprWrYOnn5a932bNpN5Ax45SZ17lWMkJwllvOXrsWMkubNFCl3dysAeTn0wTIo4WZccWX3pdi8HFfz5G47GY1pjbDzBwja3DyWW7YdZc+OEHOcvboYM0UfD317O9KnEREdL9qnBhKexTv/7ds70ZVMdZV1+yl6wXhB991NEjUJlAXFKUaULkqUKEbPEl+kJenBsuoVi7d7gQfQwnuw+xZjRgw8lusHDVTzxz8HP5Ah07Vr44vb0d/VZUVnDsmCw3f/IJXL8u1dDatQNPTylTmkG0BWD2k/WCsFJIUtTI2Wc5v+86Me4bMCqtx6nXEmJyHaWwpx+e196m+MmCND61DGvs37Q5Hou7LT9vPdOfN5dO1CItKuleflkym52cZKn5+edl6dkB0qoFoM6mMw8NwirL2bMHFr/pxendZ6BPR7BGYRqQz60Yc9ouo9sJd6589gFF/t5KlNWJHys1ZVqLdhwpVUWOF2kAVg9z6ZLMeAcPluSqRo1kxjtgAJQo4dChpUULQJ1NZy4ahFWWceSINLj68kvwqLSNAgN6co0oACwYjIypTvenXoczZyhSsiQHnx/NKM96HIx1o4SnO1P1al8lxDRh2zY5z/bNNxAdLf3Eu3aVxi6ZRFqcTU+r2bRKG8kOwoZh5Ae+AIoAO03THJzmo1LZXnKWw86ehbfekiZFLt57qfDmeI5afsDT1RPnaCfsdhsusSatFv0CFVrI0uFTT1HFyYk1Gfy+VBYUHi7Ly//8I0fUBg+Wus5Vqjh6ZP+RFmfT02I2rdJOSmbCvYAvTNP8wjCM5YZh1DVNMw3PIKnsLqnLYZcuwdSpUvHPzH+UcqMncNTlSy655mWyxzO8+NVp9p36m00VnfGv1o6GP7+tdcFV0hw5IjPf3r0hd24pbztoUKY/opYWLQC10lvmkuxzwoZh9ACqAdOBH4AupmleTOjxaX5OWGV5CZ3x9fJ0548xLQkJgXfflVbP4T6rKNRhCtdc/8HNyZUXo2sz6tMj5D97Wc6KP/+8lAL09HTAO1FZis0GP/4oS86//irlbc+fl9lvDpJYoRuVeul9Tngr8CTwAnAIuJaC51A5WELLXkGXopk2DaZPh5CYS5QcPILwvF9zBXAyDb5aFkW7w3/Ck0/Cx8OhdWs9K66SZtMmuVg7cwa8vKSV6YABOS4AQ9rMplXaSUkQngAMMU3zhmEYLwN9gYX3PsAwjEHAIAAfH59UD1JlPQ/b831wOcyMtXBzjzdh2yowNjaE8r3fJcp7FsG2cAwTTANM02Rfx8a0G7BEEmaUepi4RCsXF6hTRz4zFSveyRfAKWfnpHao5aVBN5NIyTQiP+BnGIYVqA//7eZmmuZC0zTrmqZZt3Dhwqkdo8pi4pa7gkMiMLm757tqVzAgySXuzlZMu0HYvpIEf9yc65vLkP/JGeR5vRTHSkzhqYPRLP8G3GwGViy4uLjjP3SGBmD1cLduSSnSOnXkaNGUKXK7j48sQXfqlOMDsMpcUhKEpyIz31CgALAiTUeksryHHYEAeKqGF+3y1OPCT25cDfs/jEdfJd9YLy5VfpPmh26w62MrX9o7UeLllZR3+5A8MT0pb0zn4hVdVVEPMWOGLDUPHAixsVLh6vPPHT0qpR4q2ZeEpmn+DWgKqkpQQnu+wf/f3p3H2Vz2fxx/XbMa2yAkY8hSIu7KPiFLirRYbkUiLdpIm5+6RbRJtGdQpFRK6a50504qpSSpQSWJ4rbvDGMZs53r98c1wwwznBkzc5Z5Px8PDzNzvnPm+p7vme9nru3zSUxm7lwYMQKW7/kDM6ALhKSRbqDOVpi4pBJx3YfAs7cze4fNXDxSh2jqkJSEEgpIThkZ8Pnnbm1ARASEhkLnzm6xXps2ygUuAUGrWqTQ5bbV4cimSuyd1ZquV3rIKPssFfpd5oorGAix0KvxdcT9tA0efRSqVz9lb1pKsF27XIHounVd/ubZs93Xhw51mVzatlUAloChICyFLmvOFyBle3l2zGrOrneb0Tz6aRreU57fLh1GpYMpRNgQQgklMjyKDt3uy1H+TQkF5AQHD7oVzrGxrgBH7douu1WPHr5umUiBaYWCFLruF8WwaV0YI15OwOP5nHZVn2Zvp6/4OiaVeklhvBNxPb1HvshPR9ayYP0C2p/dnrjYuBzPoYQCArjSgStWQIsWLqnGmjVw662ubq8Ss0gQUBCWQvW//8Fjj0HC15+SPOAO0kM9zDNQOTWcKbGDuanfs4RHlAIgjqonBN8shZGeTwLY//4HkyfDtGlukdWWLS6T1Q8/aKhZgoqCsBSKbdvg6UePcGDaLLpUHs9P/1xJeua7KwTDPZc/wm3tHvH6+ZRQoIT69Ve3cu+zz1wilu7d3UKrMmXc4wrAEmQUhOW07NkDU0ZuJOy1V+hb/hUmdEukb2OIMhGEhXiw1hIRGkGnOp3y/dxKKFBCJCa6/b0xMW7Fc0ICjBzpcjnXqOHr1okUKQVhKZADSZaP711AxRkT6Ft6NmOusDx8kSE8LJKH4u5jWOsHWb17dZ5zviL88ovL4/zOO9C7tyuT1aSJK5sVHu7r1okUCwVhOcHJUk4e2X2QRXe+TczseCrX/oMnBoSxtIaB0DDuanYnD7d9mGplqwEQFxun4Csnmj3bVehYtAiioqBfPzfknEUBWEoQBWHJIa8yg6XWrqPG5H9T85vpNItI4qZulZjdGCAdTAgPNpnIuCtu82nbxY9t2QLVq7s53W+/hR074Pnn4aaboGJFX7dOxGcUhCWH7EkyQjwZtFu7lP4Lv6TjrsUkhodx/2WNeC/ubw6bvS5ruAEsvLHkB+Kqd9UcrhxjrateNHGi6/3OmweXXgpPPgnPPacKWCIoCMtxtu5LJjr5ANf+9iX9f/qcWoe3si70LLq2vpwlXX5hb8YvVAhpSVRyHHsjJmFtOoYwQtPP55l5qxWEBVJS3NaiiRPhjz+gUiWXzercc93jWSudRURBWLL59Vde+noSnZbN59ezUniyaS0WRV3N6qY/44n8gktqXMJTHT+i36R9WCA8NYYjISso5WlMpKeBslmVdElJrj6vMa5eb0yMW2zVu7eb+xWREygIB6CTLZzKt7Q0+OgjiI+H77/nspAoRsdczrMDPseGbgCzgQhbg+Gt3mH05ddjjKF6ha/Zsi+ZSE8DIj0Njj6VslmVQOnpMGeOe/+sWQPr1rn0o8uXw5lnal+vyCloUibAnKpWr9e2b4fHH4dataBPH7Yv28r9PEvdC6cz8eafsWGuuAIYrm1wA4927ovJvKFmzw2dRdmsSpjdu48VUejRwwXgO+90f9QBVKumACziBfWEA8zJqgudsjdsLSxezKYnxlPty/8SnpHOV+U68gJT+a5WGSr1Hsk2FlG9XHXSD0eQ4ckgIjSCwRd3y/E0ymZVgqWluS1EP//siihceim8+CJcfTWE6XYikl/6rQkwBaoulJwMM2e6IcPly4mOKMO/GnRhesVz2Hu4IWFx40ivvJDy5aoz+ZLJ3HLRLSzduvSkiTaUzaoEOXIE3n/fLbTq0AHGjXN1e1etgvPO83XrRAKagnCAyVd1oexJ8PfuJb1BI4bXG8uE1GakdLsawuaAgQxbmlqhd7BqyAtEhbvnUaINYcMGeOUVeO01N/x83nnQsKF7LCREAVikEGhOOMCccj7W44Evv4Ru3dx83fPPk9a2I6/fuICKG3/j2T3XktbjPgg7krnH11A+vRvm4NVHA7CUYNYe+/jhh2H8eGjbFubPd9uNBgzwXdtEgpB6wn7mVCuf85yPrVsOJkxwQ4arV0OVKqQNe5jXwu5k5Cs12Ju2hTq3DuJgxal4CAFCwVoMYUR5mmhlc0mXlARvvQWTJrmh58aN3TajsWOhZk1ft04kaCkI+5G8UkYCJwTio5//8QdMfMrdQA8ehJYtSX/jbV7bfy2Pj4tk2/7d1O7/fxyKmcgmMuhSsx9//X0ZhzK2Hd3jWyG0kVY2l1R//unWCrz5pnv/tGjh/geoU8e3bRMpARSE/YjXK5+z9mZOmABff+32ZV5/PRl3DuadNc0ZPRrWh3xJ9LVPE1F5MRtsCv0b92d0u9HUrlg7W2+7gVY2l2SHD0Pz5pCaCn36uCIKLVr4ulUiJYqCsB855crn3bvdIpnJk2HjRoiNhaeewt46kI+/r8LIW2DV34epfP0wqDOJ/UCoDWVGjxlc3/j6o8+nlc0l1K5dbpHeDz/AJ59A6dLwwQfQtClUqeLr1omUSFqY5UfympftcGiTqzZTo4bbm1mvHnz4IXbtOuY1GU7zrlX453Wp7K07iQqj67K7zqQc379+3/qib7z4r4QE9/6JjXXvn0OHYP9+91iXLgrAIj6kIOxHsq98jkhPo9vKb5g94/94Pf4u+Pe/4ZZb4PffYf58Fp3Zk/adwuhyRQbrK7xJ5cfqs6PZYBpVr8ekrpOICosi1IQSERpB+7Pb+/bExHc++cQNOWe9f1audCudK1TwdctEBA1H+5XuF8VQauc2to1/iat/nEPlw/s4WLM2vPSS2xoSHc3y5TDySvhsxSIi206i0mOL2OPZQJMqTXi742Q61+3MJ79spZ5JZ1NqArGlm7Fjd02I9fXZSbHYvBlefdWlIx040CXViI+Hfv0gOtrXrROR4ygI+wNrYeFCiI+ny0cfub2+XbvCkCGUvewyCAlh9WoYdTvMmmWJ6vAC3Pp/pGBJ9RjGdBzD8DbDMcZkW2Fdh2jqkJREriusJYhYC99954Ltxx+798+gQe6xUqXcgisR8UsajvalQ4dgyhS44AJo1w6++gruvx/+/tutfu7cmQ2bQrjlFpeo6D+/LCR29CUktxsKuKQKISYEgzlaXOFkK6wlSN1xB7Rv74aZH3gA1q51AVlE/J6CsC+sXeuKnNeo4W6gISFu1fPmzfDMM1CnDtu3wz33uDroM75eSs3hV3Ck7yWkl1vL0Lihec75Fii3tASWrPfPlszKWX37Hnv/jB8PtWv7tn0i4jUNRxcXjwfmzXM9lLlzITQUevWCu++Giy8+WvYtMRHufXYx7y1eQPqeWpx9/0f8L+pDkqIqMb71eAa3GEzp8NL8s8E/cy2wkK/c0hI4cnv/tGwJ113nesHt2/u6hSJSAArCRW3fPnjjDZcO8O+/XZ3VUaPg9tuhevWjhx08CC+/DE+9vZhD/+wIl7jcztvDohh18SgeiHuA6FLHFtbkVWBhWOf6ObJugWr9BryUFLjwQpfdKo/3j4gEJgXhorJiheu1zJjhMhO1bu1y8fbs6TJcZUpJcYtZx4yBnclbqXDH/RB+BACD4YG4B3isw2Ne/1jV+g0SWVuJ7rkHIiPh2mtd1aJevXK8f0QksBmbvWpKEWjWrJlNSEgo0p/hN9LS3L7M+Hj49lu3MvWGG9zq1IsuynFoerpL9/zoo7Bp9x5i+45jR60JZNi0o8dEhEYw/8b5KilYUqSnw6efunSk33zj3j8bNyqZhkiAMcYstdY28+ZY9YQLw44dMHWqq726ZQtbK1bjzfY3812bq7mjZ4scvdCPlm5h+HOJrJt3NumHPFS65llK/+NFNmcc5IZGN/Bou0fZeWhnrvO9EsSWLHHzuxs3uqpFTz8Nt94KlSv7umUiUoQUhAvKWvjpJ9drmTUL0tLY0aodj7YeyLxaTfCEhELasT263S6M4ZEJu3n2qXKklPoLc81QqLGIvWGHaFmlC6/1eIZGVRsBULdSXQXfkmDZMjhyxC3Mq1fP7UN76SW46ioI06+mSEmg3/T8OnLE1VuNj3c5ecuVg7vugkGD6PnxlhNWJienZTBy8g6eWRnDD0vKY7o8CM1ewhrAGiql3ENoYrejAViCXGoqfPihe//88AN06OAqYZ1xhlv1LCIlSoGDsDFmEjDXWvtpIbbHf23c6Iabp0511YwaNICJE6F/fxeIga37/s7xLSnbotn3XX02bKhExUveJmTYw3hKbc7KswEYPCH7tIe3pJg61a1s3r7d9XxffNEVVhCREqtAQdgY0xaoFvQB2Fq3QCY+3i24AujWzS206tjx6N7eLFl7dFN3l2X/wnM5vKYa5oIPCB82gsTSf1PG1CMi5Rb2hc/A2nQMYZTyNNYe3mBlLfz4IzRq5P5Qs9Yt0BsyxOV0DlGuHJGSLt9B2BgTDkwFPjPGdLPWflL4zfKxgwfd0uX4eFi1yi2OeeghuPNOt2gmD/3Pb8hDIzJISlwPzYcTcs1CPGXXUbVsPV7oMovQlFaM+HglkakNOBKyglKexlQIbaQ9vMHmyBF47z23XmDZMjdiMmgQ3Hab298rIpKpID3hG4E/gPHAEGNMTWvthOwHGGNuB24HqHmSoOV3Vq92STWmT4ekJFfsfPp06N3bbRfJw9at8OSTMHVqNWyTV6DnYDAePECvc+5hZp/nCAtxL3WICeGZeRFs3ddAe3iDTUYGPPLIsSmLhg3d+6l/f/f4cSMnIiIFCcIXAVOstduNMTOAMUCOIGytnQJMAbdP+LRbWZQyMuCzz1yv5csvITzcBd2774YWLU5649y9G8aNcx3mtDOWEzNsJBsjPzv6eKgJpUlstaMBGFwyDQXdIGItrFkD9eu7VJKLFrnELEOG5DplISKSXUGC8N9AncyPmwEbCq85xWjPHnj9dddTWb8eYmJcd3bgQDjzzJN+a1ISvPACPPccHIj8k1p3j2JD2Q84UKoigxoN4o1f3iA1I/WE4goSRA4fhnfecX+8rV4NmzZB1aquElZ4uK9bJyIBoiBBeBrwujGmDxAO9CrcJhWxX35xXdd33nFzd+3aucpF3bqd8uaZnAwPTljM6/MXcHhjfWrePIdDld5kT0RpHmn1CEPjhhJdKpp+/+inZBvBavt2eP55V7UoMdGVoZw06egKeQVgEcmPkpG2MjUVPvrIBd9Fi6B0aTdPN3gwNG58ym9PS4Np0+CRVxez+8qOEJYCxhIeEs6QFkP4V5t/UaWMUgsGLWth/36oUMGVEWzQALp3d0PObdpoyFlEclDayixbt8KUKa5CwvbtULeu68XcdBNUrHjKb8/IgJkzYfRoWLdtL+Vv+9cJxRWe7vR0EZ+E+MzBg/D2227I+Zxz3Da1unVh2zaXXENE5DQFXxC21vV24+NdZqL0dOja1S208nJvprXufjtyJKz86wDVur9E6cbPkJSRRKgJBVxxhW71uxX12YgvrF3r3j9vvOF6wE2buupFWRSARaSQBE8QPnzYdVvj4928b3S0Gy4cNMhlJ/KCta563IgR8NOyI1S5YjLl+jzF9ozddD+nO090eIIDKQc03xuMPB73f0iIWy8QH+/KBw4ZAq1aachZRIpE4M8Jr1sHkye7SdvERJedaMgQV0KwTBmvn2bxYrj9nhR+3/kTtH2GkLO/xxOeyGV1LuPJjk/SIqZF0Z2D+M6BA/Dmm27IecwY1+NNTHSL9s46y9etE5EAFPxzwh6P2woSHw9z5rjeS8+ebsi5bdt89Vp++80NO386xwOXPglXjgFj8VhD1Yy7GNRoBC1itK836KxZ494/06e7QNyypVt4BV6tFxARKQyBFYT373e9lokT3U20alU3dnzHHVCjRr6e6q+/3IKrmTMtZZr8h/BhD5JWek2O4gqpnkM8M2+1kmsEG2vhmmvcKErv3m7kpIVGOkSk+AVWEH7rLbj3XjdHN2OGGzqMjMzXU2zaBE88AdNet4SfO5/qo0awNeQnwjwxRKf1IylsVo7iCqpwFAT273c93nffdQU5Spd276VatU6ZmEVEpCgFVhAeMADi4qCZV0PtOezcCWPHQvwni/E0eoszhv/ErrBlhJaPZVr7abz2eSxb96cS5bngaHGFSE+DHBWOZi/fwjPzVrN1X7LyPgeCP/90Q85vvum2G8XFue1Fdeuq5ysifiGwgnD58vkOwPv2ufSSL7wAh855C3PjzVjjYRdwf8v7GdtpLJFhkVRiC8M/WgFpDYj0NAAgKjz0aIWj2cvd48lpGQBs2ZfsjgcFYn/0xx9w/vkQEQF9+rgh5wL88SYiUpQCKwjnw6FDrhM0bhwkhqwh9vZRHIp+/+iUb6gJpUqZKkSGueHsrECaV0/3mXmrjwbgLMlpGZoz9hf797t9vUlJMGqUy2r16qsus1XVqr5unYhIroIuCKemukpyTz4J25M3EtvvcZKqTmdvWCkGNBzArJWz8iyucLIKR3nNDWvO2MdWrTo25HzoEFx2mVt4ZYxq94qI3wuaIJye7tZqDZ+0mO3l5lC5x1+EV/uEHSEwpNkQhrcdTtUyVbmj6R0FSrZRvUIUW3IJuNnnjKWYvfyyW6gXGQnXX++GnJs08XWrRES8FvBB2ONxtRkeeQT+TJ0H/a6CkHR2A1efczXxXeOpGV0TyFpYlczWff9gToVkhnXe4vVQ8rDO9XPMCUPOOWMpBvv2uSHn1q3dwqrOnV2CjdtugyoqoCEigefUiZT9lLXw+efQvDlce8NB9jR8ivAB3SAkHXBzvnE14nIE4OEfrWDLvmQsxxZWzV6+xauf1/2iGMb2bExMhSgMEFMhirE9G2s+uDisWuUqXtWoAQ88AJ9+6r5evz48/LACsIgErIDsCS9c6HJ0LFx8hEqXv0r57k+xy7OTNrFt+HnLz6R70k+Y8y2MhVUnmzOWInLDDW5/b0QE9O2rIWcRCSoBFYR/+w3+9S+YOy+d6HbTqTDqMfZ6NtOxVkfGdBxDqxqtWLxpca5zvlpYFSD274f33oOBAyE01A11NGzohpy1yllEgkxABeH/LF/EVxUnED1qEfvZTMuzWjKm43QurXMpcPI5Xy2s8nNZiTWmT3ernM87D9q1g/vu83XLRESKTMDMCS/etJjHN3Yg7dz3SWIL4zuNZ/Gti3ME4JPN+Q7rXJ+o8NAcz6mFVX5g507o0sXt65061aUiTUhwAVhEJMgFTBBesH4BHutqvoaYENI96Zhs1ZJONucLWljlV5KS4Mcf3cdnnOFSSj7xhEvsPX06NG3q0+aJiBSXgBmObn92eyJCI/JMtOHNnK8WVvlYVvnAN96AqCjYvNktuPr+e1+3TETEJwImCMfFxjH/xvl5JtrQnK8fW7rUbeSeOxfCw10u53vucQFYRKQEC5ggDC4Q55XlSsk0/MyBA5CSApUru+Hm5cvhscdcKslq1XzdOhERvxBQQfhkTlWAQYrJ2rVuyPn11+HGG2HCBLjkEtiwQT1fEZHjBE0QBs35+tQ338Dzz8N//+v29153HfTv7x4zRgFYRCQXAROE3R5g9XL9SnKyW2AF8M47sGQJjBwJd94J1av7tm0iIgEgILYonW7eZylk69fDsGEu0C5Z4r729NOwcSM8/rgCsIiIlwIiCJ9qD7AUA2thwQLo2RPq1oUXXoDLL4dy5dzjlStDqVI+baKISKAJiOFo5X32IWvdnG5yMvTo4eZ7H3oIBg1yVY1ERKTAAiIIaw+wD2zeDJMmuZJV334LpUvDF19Ao0bH5oFFROS0BMRwtPI+FxNrYdEi6N0bzj4bxo1zw8z797vHmzdXABYRKUQB0RPWHuBiMncuXHklREe76kWDB0Pt2r5ulYhI0DLW2iL9Ac2aNbMJCQlF+jOkgLZvh1degSpVXMBNS4O33nI94bJlfd06EZGAZIxZaq1t5s2xAdETlkKWkAAvvQTvv+8C7623uq+Hhx/7WEREilxAzAlLIXrwQTe3O3u2S6qxZg289pqvWyUiUiIVuCdsjDkT+Nxae1EhtkcK2+7dMHWqq1xUuzZcc41LpnHzzW7uV0REfOZ0hqOfBbRU1l+tWAEvvwwzZsCRIy7gDhoEbdq4fyIi4nMFCsLGmI7AIWB74TZHTpvHA1dd5VY6R0W5SkZDhrj9vSIi4lfyHYSNMRHAI0APYHYex9wO3A5Qs2bN02mfeGP/fhd0+/SBkBBo3BjatYOBA+GMM3zdOhERyUO+tygZY0YBq6y1HxhjFlhr25/seG1RKkJ//eWGnKdPh4MH4c8/ob4SmIiI+FJ+tigVZHV0J2CwMWYBcKExRktri9v69S6pxrnnwpQprqhCQoICsIhI+naGxgAACSJJREFUgMn3cLS19pKsjzN7wgMLt0mSq0OHYNMmOO88qFTJ9YJHj3bbjKpV83XrRESkAE4rWcephqKlEGzcCBMnum1GNWrAr79C+fKwerWrbiQiIgFLGbP81bJlMHYsfPyxK6zQsyfce++xxxWARUQCnoKwP0lJgYwMVzZw5UqYPx+GDnV5nbXKXEQk6ChtpT/YsQMefxxq1YL4ePe13r3dHPC4cQrAIiJBSj1hX1q+3BVSmDkTUlPhiiugVSv3WESE+yciIkFLQbi4WXtsPvehh+CHH+C221xWK20xEhEpURSEi8u+fa5a0auvwtdfQ2wsTJ7sMlpVqODr1omIiA9oTriorV7tFlbVqAHDhkFMjAvIAHXrKgCLiJRg6gkXpT17XB5nY6BvX7fF6MILfd0qERHxEwrChenwYXj7bZdQY9IkN9T87rvQti2ceaavWyciIn5GQbgwbNrkslpNmQKJidC0qQvIpUtDr16+bp2IiPgpBeHT9Z//uGxW1kKPHnDffdC6tTJaiYjIKSkI51daGvz731CxInTp4oaahw6FQYNcsg0REREvaXW0t/bscbmca9d2i6ymTXNfr1jRZbVSABYRkXxSEPbG2LFui9HDD0PDhjBnDrz/vq9bJSIiAU5BODceD/z3v5CU5D6PjYUbb4Tff4cvvoArr4QQvXQiInJ6FEmyO3jQrXI+7zy46iqYMcN9vV8/l+nq/PN92z4REQkqCsLgygcOG+aGnO++283zzpzpcjqLiIgUkZK7OtpaWLsW6tWD0FBX0ahLF7fFKKuSkYiISBEqeUE4NRU++ABefNFlttq4EapVg3nzXDAWEREpJiVnODoxEcaMgbPPdnO8Bw7Ayy9D+fLucQVgEREpZsHfE05JgchIV7lo9Gjo1Mnt8e3cWSucRUTEp4IzCHs88Pnnbsg5IsLt661dGzZscKUERURE/EBwdQUPHnTVixo0cHt5V66ENm3cIixQABYREb8SXD3h+HgYPhyaN3clBHv1gvBwX7dKREQkV4HdE168GHr3dqudwe3rXbQIliyB669XABYREb8WeEE4LQ3ee8/t5b34Yre1aPdu99gZZ7ivqYygiIgEgMAbju7aFb76Cs45xw0/DxgAZcv6ulUiIiL5FnhB+L773L8rrtAWIxERCWiBF4SvvNLXLRARESkU6kqKiIj4iIKwiIiIjygIi4iI+IiCsIiIiI8oCIuIiPhIvldHG2OigfeAUOAQ0Ntam1rYDRMREQl2BekJ3wA8b629HNgOdCncJomIiJQM+e4JW2snZfu0CrCz8JojIiJSchR4TtgYEwdUtNb+WIjtERERKTEKFISNMZWACcAteTx+uzEmwRiTsGvXrtNpn4iISNDKdxA2xkQAHwDDrbUbcjvGWjvFWtvMWtusSpUqp9tGERGRoFSQnvCtQBNghDFmgTGmdyG3SUREpEQoyMKsycDkImiLiIhIiWKstUX7A4zZBeQ6bF1AlYHdhfh8vqRz8U/Bci7Bch6gc/FHwXIeUPjnUsta69VcbJEH4cJmjEmw1jbzdTsKg87FPwXLuQTLeYDOxR8Fy3mAb89FaStFRER8REFYRETERwIxCE/xdQMKkc7FPwXLuQTLeYDOxR8Fy3mAD88l4OaERUREgkUg9oRFRESCgt8GYWPMNGPMYmPMyNM5xteMMdHGmLnGmC+MMR9nZhw7/pgwY8zGzOQnC4wxjX3R1pPxto3GmMeMMT8bYyYWdxu9ZYy5K9t5/GKMeTWXYwLhmpxpjFmY+XG4MeZTY8wiY0yu6WTzc1xxO+5cama+5l8bY6YYY0we3xNjjNmc7Rr5RXq+487F6zb62/3suPN4LNs5/GmMGZ7H9/jdNcntHuzta10c18Qvg7AxpicQaq2NA+oYY84pyDF+wpvSj/8AZlpr22f+W1GsLfTOKdtojGkKtAFaADuNMZ2Ku5HesNZOzjoPYCEwNZfD/PqaGGMqAm8CZTK/NARYaq1tDfQyxpTL41u9Pa7Y5HIudwB3WWs7ArFAXn8AtQTGZLtGPk9Un8u5eNVGf7ufHX8e1trR2X5nfgfeyuNb/e6acOI9uA9evNbFdU38MggD7YFZmR9/gbuxF+QYn7PWTrLWfpn5aV6lH1sBVxljfsr8yyvfmcyKgTdtbAd8aN1Cg3lA22JtYT4ZY2KAM621Cbk87O/XJAPoDSRlft6eY78P3wF57Xn09rjilONcrLUjrLWrMh87g7yTKLQCBhpjlhljnir6Znrl+OvibRvb41/3s+PPAwBjTHNgs7V2Sx7f53fXJJd7cD+8e63be3ncafHXIFwGyLrIe4EzC3iM3zAnL/34M9DJWtsCCAe6FmvjvONNGwPqmgCDyTsFq19fE2ttkrV2f7Yvefva+901yuVcADAuL/1Ka+3WPL51Lu5G2RyIM8b8o+ha6Z1czsXbNvrVdcnrmgD34iro5cXvrkmWrHswsAk/+l3x1yB8EIjK/LgsubfTm2P8gjlF6UfgN2vttsyPEwB/HFr3po2BdE1CgA7AgjwOCYRrkp23r31AXCNjTB3g/4D7TnLYD9baA9baDGA5/nmNvG2j318XY0wFoKq1du1JDvPLa3LcPdivflf87kJnWsqxrv8FwPoCHuNzxovSj8DbxpgLjDGhQHfg12JroPe8aWNAXJNMbYElNu89eoFwTbLz9rX3+2uUOR85E7glj95YlnnGmLOMMaWBy3Fzlf7G2zb6/XUBugGfneIYv7smudyD/et3xVrrd/+A8rib3vPAqswX4MlTHBPt63bncS53AYm4HtcCYHQu59II+A1YgVvU4PN253IeOdoIVAJeO+6YEGAR8BKwGqjt63af5HyeAnpmftwwEK9JZjsXZP5fC1iZ+dr/DIQCHYG7jzv+hON8fQ65nMs4YFu235l2eZxLB+DPzOt0d3G318tzOaGNebzf/PJ+lnUemR+/CzTJ9nlAXJNc7sEDjn+tfXlN/DZZR+Zfw5cB31lrtxf0GClexpgo4EpgmbV2na/bU5IYY6rj/nKfZ0/Sg/T2OCleup8VH29f6+K4Jn4bhEVERIKdv84Ji4iIBD0FYRERER9REBYREfERBWEREREfURAWERHxEQVhERERH/l/oC4UecuEqXoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(x1, y2, 'o', label=\"data\")\n",
    "ax.plot(x1, y_true2, 'b-', label=\"True\")\n",
    "ax.plot(x1, res2.fittedvalues, 'r-', label=\"OLS\")\n",
    "ax.plot(x1, iv_u, 'r--')\n",
    "ax.plot(x1, iv_l, 'r--')\n",
    "ax.plot(x1, resrlm2.fittedvalues, 'g.-', label=\"RLM\")\n",
    "legend = ax.legend(loc=\"best\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
